26 #include <queso/VectorSequence.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
33 template <
class V,
class M>
36 unsigned int subSequenceSize,
37 const std::string& name)
39 m_env (vectorSpace.env()),
40 m_vectorSpace (vectorSpace),
42 m_fftObj (new
Fft<double>(m_env)),
44 m_unifiedMinPlain (NULL),
46 m_unifiedMaxPlain (NULL),
47 m_subMeanPlain (NULL),
48 m_unifiedMeanPlain (NULL),
49 m_subMedianPlain (NULL),
50 m_unifiedMedianPlain (NULL),
51 m_subSampleVariancePlain (NULL),
52 m_unifiedSampleVariancePlain(NULL),
54 m_unifiedBoxPlain (NULL)
56 if (subSequenceSize) {};
59 template <
class V,
class M>
63 this->deleteStoredVectors();
64 if (m_fftObj != NULL)
delete m_fftObj;
68 template <
class V,
class M>
72 unsigned int unifiedNumSamples = 0;
74 bool useOnlyInter0Comm = (m_vectorSpace.numOfProcsForStorage() == 1);
76 if (useOnlyInter0Comm) {
77 if (m_env.inter0Rank() >= 0) {
78 unsigned int subNumSamples = this->subSequenceSize();
79 m_env.inter0Comm().template Allreduce<unsigned int>(&subNumSamples, &unifiedNumSamples, (int) 1,
RawValue_MPI_SUM,
80 "BaseVectorSequence<V,M>::unifiedSequenceSize()",
81 "failed MPI.Allreduce() for unifiedSequenceSize()");
85 unifiedNumSamples = this->subSequenceSize();
92 return unifiedNumSamples;
95 template <
class V,
class M>
99 return m_vectorSpace.dimLocal();
102 template <
class V,
class M>
106 return m_vectorSpace.dimGlobal();
109 template <
class V,
class M>
113 return m_vectorSpace;
116 template <
class V,
class M>
123 template <
class V,
class M>
131 template <
class V,
class M>
135 unsigned int numPos = this->subSequenceSize();
137 this->resetValues(0,numPos);
138 this->resizeSequence(0);
144 template <
class V,
class M>
148 if (m_subMinPlain == NULL) {
149 m_subMinPlain = m_vectorSpace.newVector();
150 if (m_subMaxPlain == NULL) m_subMaxPlain = m_vectorSpace.newVector();
151 subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
154 return *m_subMinPlain;
157 template <
class V,
class M>
161 if (m_unifiedMinPlain == NULL) {
162 m_unifiedMinPlain = m_vectorSpace.newVector();
163 if (m_unifiedMaxPlain == NULL) m_unifiedMaxPlain = m_vectorSpace.newVector();
164 unifiedMinMaxExtra(0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
167 return *m_unifiedMinPlain;
170 template <
class V,
class M>
174 if (m_subMaxPlain == NULL) {
175 if (m_subMinPlain == NULL) m_subMinPlain = m_vectorSpace.newVector();
176 m_subMaxPlain = m_vectorSpace.newVector();
177 subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
180 return *m_subMaxPlain;
183 template <
class V,
class M>
187 if (m_unifiedMaxPlain == NULL) {
188 if (m_unifiedMinPlain == NULL) m_unifiedMinPlain = m_vectorSpace.newVector();
189 m_unifiedMaxPlain = m_vectorSpace.newVector();
190 unifiedMinMaxExtra(0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
193 return *m_unifiedMaxPlain;
196 template <
class V,
class M>
200 if (m_subMeanPlain == NULL) {
201 m_subMeanPlain = m_vectorSpace.newVector();
202 subMeanExtra(0,subSequenceSize(),*m_subMeanPlain);
205 return *m_subMeanPlain;
208 template <
class V,
class M>
212 if (m_unifiedMeanPlain == NULL) {
213 m_unifiedMeanPlain = m_vectorSpace.newVector();
214 unifiedMeanExtra(0,subSequenceSize(),*m_unifiedMeanPlain);
217 return *m_unifiedMeanPlain;
220 template <
class V,
class M>
224 if (m_subMedianPlain == NULL) {
225 m_subMedianPlain = m_vectorSpace.newVector();
226 subMedianExtra(0, subSequenceSize(), *m_subMedianPlain);
229 return *m_subMedianPlain;
232 template <
class V,
class M>
236 if (m_unifiedMedianPlain == NULL) {
237 m_unifiedMedianPlain = m_vectorSpace.newVector();
238 unifiedMedianExtra(0, subSequenceSize(), *m_unifiedMedianPlain);
241 return *m_unifiedMedianPlain;
244 template <
class V,
class M>
248 if (m_subSampleVariancePlain == NULL) {
249 m_subSampleVariancePlain = m_vectorSpace.newVector();
250 subSampleVarianceExtra(0,subSequenceSize(),subMeanPlain(),*m_subSampleVariancePlain);
253 return *m_subSampleVariancePlain;
256 template <
class V,
class M>
260 if (m_unifiedSampleVariancePlain == NULL) {
261 m_unifiedSampleVariancePlain = m_vectorSpace.newVector();
262 unifiedSampleVarianceExtra(0,subSequenceSize(),unifiedMeanPlain(),*m_unifiedSampleVariancePlain);
265 return *m_unifiedSampleVariancePlain;
268 template <
class V,
class M>
272 if (m_subBoxPlain == NULL) {
276 this->subMaxPlain());
279 return *m_subBoxPlain;
282 template <
class V,
class M>
286 if (m_unifiedBoxPlain == NULL) {
289 this->unifiedMinPlain(),
290 this->unifiedMaxPlain());
293 return *m_unifiedBoxPlain;
296 template <
class V,
class M>
301 delete m_subMinPlain;
302 m_subMinPlain = NULL;
304 if (m_unifiedMinPlain) {
305 delete m_unifiedMinPlain;
306 m_unifiedMinPlain = NULL;
309 delete m_subMaxPlain;
310 m_subMaxPlain = NULL;
312 if (m_unifiedMaxPlain) {
313 delete m_unifiedMaxPlain;
314 m_unifiedMaxPlain = NULL;
316 if (m_subMeanPlain) {
317 delete m_subMeanPlain;
318 m_subMeanPlain = NULL;
320 if (m_unifiedMeanPlain) {
321 delete m_unifiedMeanPlain;
322 m_unifiedMeanPlain = NULL;
324 if (m_subMedianPlain) {
325 delete m_subMedianPlain;
326 m_subMedianPlain = NULL;
328 if (m_unifiedMedianPlain) {
329 delete m_unifiedMedianPlain;
330 m_unifiedMedianPlain = NULL;
332 if (m_subSampleVariancePlain) {
333 delete m_subSampleVariancePlain;
334 m_subSampleVariancePlain = NULL;
336 if (m_unifiedSampleVariancePlain) {
337 delete m_unifiedSampleVariancePlain;
338 m_unifiedSampleVariancePlain = NULL;
341 delete m_subBoxPlain;
342 m_subBoxPlain = NULL;
344 if (m_unifiedBoxPlain) {
345 delete m_unifiedBoxPlain;
346 m_unifiedBoxPlain = NULL;
352 template <
class V,
class M>
356 unsigned int initialPos,
363 this->deleteStoredVectors();
364 unsigned int currentSize = this->subSequenceSize();
365 this->resizeSequence(currentSize+numPos);
367 for (
unsigned int i = 0; i < numPos; ++i) {
369 this->setPositionValues(currentSize+i,tmpVec);
375 template <
class V,
class M>
381 if (m_env.subDisplayFile()) {
382 *m_env.subDisplayFile() <<
"Entering BaseVectorSequence<V,M>::subPositionsOfMaximum()"
383 <<
": subCorrespondingScalarValues,subSequenceSize() = " << subCorrespondingScalarValues.
subSequenceSize()
384 <<
", this->subSequenceSize = " << this->subSequenceSize()
390 double subMaxValue = subCorrespondingScalarValues.
subMaxPlain();
393 unsigned int subNumPos = 0;
394 for (
unsigned int i = 0; i < iMax; ++i) {
395 if (subCorrespondingScalarValues[i] == subMaxValue) {
400 V tmpVec(this->vectorSpace().zeroVector());
403 for (
unsigned int i = 0; i < iMax; ++i) {
404 if (subCorrespondingScalarValues[i] == subMaxValue) {
405 this->getPositionValues (i,tmpVec);
411 if (m_env.subDisplayFile()) {
412 *m_env.subDisplayFile() <<
"Leaving BaseVectorSequence<V,M>::subPositionsOfMaximum()"
419 template <
class V,
class M>
425 if (m_env.subDisplayFile()) {
426 *m_env.subDisplayFile() <<
"Entering BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()"
427 <<
": subCorrespondingScalarValues,subSequenceSize() = " << subCorrespondingScalarValues.
subSequenceSize()
428 <<
", this->subSequenceSize = " << this->subSequenceSize()
435 double subMaxValue = subCorrespondingScalarValues.
subMaxPlain();
440 double unifiedMaxValue;
441 std::vector<double> sendbufPos(1,0.);
442 for (
unsigned int i = 0; i < sendbufPos.size(); ++i) {
443 sendbufPos[i] = subMaxValue;
445 m_env.inter0Comm().template Allreduce<double>(&sendbufPos[0], &unifiedMaxValue, (int) sendbufPos.size(),
RawValue_MPI_MAX,
446 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
447 "failed MPI.Allreduce() for max");
453 for (
unsigned int i = 0; i < iMax; ++i) {
454 if (subCorrespondingScalarValues[i] == unifiedMaxValue) {
461 V tmpVec(this->vectorSpace().zeroVector());
464 for (
unsigned int i = 0; i < iMax; ++i) {
466 if (subCorrespondingScalarValues[i] == unifiedMaxValue) {
467 this->getPositionValues (i,tmpVec);
474 std::vector<int> auxBuf(1,0);
475 int unifiedNumPos = 0;
476 auxBuf[0] = subNumPos;
477 m_env.inter0Comm().template Allreduce<int>(&auxBuf[0], &unifiedNumPos, (int) auxBuf.size(),
RawValue_MPI_SUM,
478 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
479 "failed MPI.Allreduce() for sum");
488 unsigned int Np = (
unsigned int) m_env.inter0Comm().NumProc();
490 std::vector<int> recvcntsPos(Np,0);
491 m_env.inter0Comm().template Gather<int>(&subNumPos, 1, &recvcntsPos[0], (int) 1, 0,
492 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
493 "failed MPI.Gatherv()");
494 if (m_env.inter0Rank() == 0) {
499 std::vector<int> displsPos(Np,0);
500 for (
unsigned int nodeId = 1; nodeId < Np; ++nodeId) {
501 displsPos[nodeId] = displsPos[nodeId-1] + recvcntsPos[nodeId-1];
503 if (m_env.inter0Rank() == 0) {
504 queso_require_equal_to_msg(unifiedNumPos, (displsPos[Np - 1] + recvcntsPos[Np - 1]),
"failed MPI.Gather() result at proc 0 (unifiedNumPos)");
515 unsigned int dimSize = m_vectorSpace.dimLocal();
516 int subNumDbs = subNumPos * dimSize;
517 std::vector<int> recvcntsDbs(Np,0);
518 m_env.inter0Comm().template Gather<int>(&subNumDbs, 1, &recvcntsDbs[0], (int) 1, 0,
519 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
520 "failed MPI.Gatherv()");
521 if (m_env.inter0Rank() == 0) {
527 std::vector<int> displsDbs(Np,0);
528 for (
unsigned int nodeId = 1; nodeId < Np; ++nodeId) {
529 displsDbs[nodeId] = displsDbs[nodeId-1] + recvcntsDbs[nodeId-1];
531 if (m_env.inter0Rank() == 0) {
532 queso_require_equal_to_msg(((
int) (unifiedNumPos*dimSize)), (displsDbs[Np - 1] + recvcntsDbs[Np - 1]),
"failed MPI.Gather() result at proc 0 (unifiedNumPos*dimSize)");
540 std::vector<double> sendbufDbs(subNumDbs,0.);
541 for (
unsigned int i = 0; i < (
unsigned int) subNumPos; ++i) {
543 for (
unsigned int j = 0; j < dimSize; ++j) {
544 sendbufDbs[i*dimSize + j] = tmpVec[j];
548 std::vector<double> recvbufDbs(unifiedNumPos * dimSize);
551 m_env.inter0Comm().template Gatherv<double>(&sendbufDbs[0], (int) subNumDbs,
552 &recvbufDbs[0], (
int *) &recvcntsDbs[0], (
int *) &displsDbs[0], 0,
553 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
554 "failed MPI.Gatherv()");
561 if (m_env.inter0Rank() == (int) 0) {
562 for (
unsigned int i = 0; i < (
unsigned int) unifiedNumPos; ++i) {
563 for (
unsigned int j = 0; j < dimSize; ++j) {
564 tmpVec[j] = recvbufDbs[i*dimSize + j];
575 if (m_env.subDisplayFile()) {
576 *m_env.subDisplayFile() <<
"Leaving BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()"
580 return unifiedMaxValue;
583 template <
class V,
class M>
587 V gaussianVector(m_vectorSpace.zeroVector());
588 for (
unsigned int j = 0; j < this->subSequenceSize(); ++j) {
589 gaussianVector.cwSetGaussian(meanVec,stdDevVec);
590 this->setPositionValues(j,gaussianVector);
593 this->deleteStoredVectors();
598 template <
class V,
class M>
602 V uniformVector(m_vectorSpace.zeroVector());
603 for (
unsigned int j = 0; j < this->subSequenceSize(); ++j) {
604 uniformVector.cwSetUniform(aVec,bVec);
605 this->setPositionValues(j,uniformVector);
608 this->deleteStoredVectors();
613 template<
class V,
class M>
616 std::ofstream* passedOfs,
617 unsigned int& initialPos,
618 unsigned int& spacing)
620 if (m_env.subDisplayFile()) {
621 *m_env.subDisplayFile() <<
"\n"
622 <<
"\n-----------------------------------------------------"
623 <<
"\n Computing filter parameters for chain '" << m_name <<
"' ..."
624 <<
"\n-----------------------------------------------------"
629 bool okSituation = ((passedOfs == NULL ) ||
630 ((passedOfs != NULL) && (m_env.subRank() >= 0)));
631 queso_require_msg(!(!okSituation),
"unexpected combination of file pointer and subRank");
636 if (m_env.subDisplayFile()) {
637 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
638 <<
"\n Finished computing filter parameters for chain '" << m_name <<
"'"
639 <<
": initialPos = " << initialPos
640 <<
", spacing = " << spacing
641 <<
"\n-----------------------------------------------------"
650 template <
class V,
class M>
659 this->deleteStoredVectors();
669 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
670 template<
class V,
class M>
673 const SequenceStatisticalOptions& statisticalOptions,
674 std::ofstream* passedOfs)
677 computeStatistics(statisticalOptions.m_ov, passedOfs);
680 template<
class V,
class M>
682 BaseVectorSequence<V,M>::computeStatistics(
683 const SsOptionsValues& statisticalOptions,
684 std::ofstream* passedOfs)
686 if (m_env.subDisplayFile()) {
687 *m_env.subDisplayFile() <<
"\n"
688 <<
"\n-----------------------------------------------------"
689 <<
"\n Computing statistics for chain '" << m_name <<
"' ..."
690 <<
"\n-----------------------------------------------------"
695 bool okSituation = ((passedOfs == NULL ) ||
696 ((passedOfs != NULL) && (m_env.subRank() >= 0)));
697 queso_require_msg(!(!okSituation),
"unexpected combination of file pointer and subRank");
700 struct timeval timevalTmp;
701 iRC = gettimeofday(&timevalTmp, NULL);
702 double tmpRunTime = 0.;
705 std::vector<unsigned int> initialPosForStatistics(statisticalOptions.initialDiscardedPortions().size(),0);
706 for (
unsigned int i = 0; i < initialPosForStatistics.size(); ++i) {
707 initialPosForStatistics[i] = (
unsigned int) (statisticalOptions.initialDiscardedPortions()[i] * (double) (this->subSequenceSize()-1));
708 if (m_env.subDisplayFile()) {
709 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeStatistics()"
710 <<
": statisticalOptions.initialDiscardedPortions()[" << i <<
"] = " << statisticalOptions.initialDiscardedPortions()[i]
711 <<
", initialPosForStatistics[" << i <<
"] = " << initialPosForStatistics[i]
715 if (m_env.subDisplayFile()) {
716 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeStatistics(): initial positions for statistics =";
717 for (
unsigned int i = 0; i < initialPosForStatistics.size(); ++i) {
718 *m_env.subDisplayFile() <<
" " << initialPosForStatistics[i];
720 *m_env.subDisplayFile() << std::endl;
726 this->computeMeanVars(statisticalOptions,
733 #ifdef UQ_CODE_HAS_MONITORS
734 if (statisticalOptions.meanMonitorPeriod() != 0) {
735 this->computeMeanEvolution(statisticalOptions,
743 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
744 if ((statisticalOptions.bmmRun() ) &&
745 (initialPosForStatistics.size() > 0) &&
746 (statisticalOptions.bmmLengths().size() > 0)) {
747 this->computeBMM(statisticalOptions,
748 initialPosForStatistics,
755 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
756 if ((statisticalOptions.fftCompute() ) &&
757 (initialPosForStatistics.size() > 0)) {
758 this->computeFFT(statisticalOptions,
759 initialPosForStatistics,
766 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
767 if ((statisticalOptions.psdCompute() ) &&
768 (initialPosForStatistics.size() > 0)) {
769 this->computePSD(statisticalOptions,
770 initialPosForStatistics,
777 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
778 if ((statisticalOptions.psdAtZeroCompute() ) &&
779 (initialPosForStatistics.size() > 0) &&
780 (statisticalOptions.psdAtZeroNumBlocks().size() > 0)) {
781 this->computePSDAtZero(statisticalOptions,
782 initialPosForStatistics,
789 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
790 if ((statisticalOptions.gewekeCompute()) &&
791 (initialPosForStatistics.size() > 0)) {
792 this->computeGeweke(statisticalOptions,
793 initialPosForStatistics,
800 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
801 if ((statisticalOptions.meanStaccCompute()) &&
802 (initialPosForStatistics.size() > 0 )) {
803 this->computeMeanStacc(statisticalOptions,
804 initialPosForStatistics,
809 std::vector<unsigned int> lagsForCorrs(statisticalOptions.autoCorrNumLags(),1);
810 for (
unsigned int i = 1; i < lagsForCorrs.size(); ++i) {
811 lagsForCorrs[i] = statisticalOptions.autoCorrSecondLag() + (i-1)*statisticalOptions.autoCorrLagSpacing();
817 if ((statisticalOptions.autoCorrComputeViaDef()) &&
818 (initialPosForStatistics.size() > 0 ) &&
819 (lagsForCorrs.size() > 0 )) {
820 this->computeAutoCorrViaDef(statisticalOptions,
821 initialPosForStatistics,
829 if ((statisticalOptions.autoCorrComputeViaFft()) &&
830 (initialPosForStatistics.size() > 0 ) &&
831 (lagsForCorrs.size() > 0 )) {
832 this->computeAutoCorrViaFFT(statisticalOptions,
833 initialPosForStatistics,
841 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
842 if ((statisticalOptions.histCompute ()) ||
843 (statisticalOptions.cdfStaccCompute()) ||
844 (statisticalOptions.kdeCompute ())) {
846 if (statisticalOptions.kdeCompute())
848 this->computeHistCdfstaccKde(statisticalOptions,
855 if ((statisticalOptions.covMatrixCompute ()) ||
856 (statisticalOptions.corrMatrixCompute())) {
857 this->computeCovCorrMatrices(statisticalOptions,
862 if (m_env.subDisplayFile()) {
863 *m_env.subDisplayFile() <<
"All statistics of chain '" << m_name <<
"'"
864 <<
" took " << tmpRunTime
869 if (m_env.subDisplayFile()) {
870 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
871 <<
"\n Finished computing statistics for chain '" << m_name <<
"'"
872 <<
"\n-----------------------------------------------------"
880 template<
class V,
class M>
882 BaseVectorSequence<V,M>::computeMeanVars(
883 const SequenceStatisticalOptions& statisticalOptions,
884 std::ofstream* passedOfs,
891 computeMeanVars(statisticalOptions.m_ov, passedOfs, subMeanPtr, subMedianPtr,
892 subSampleVarPtr, subPopulVarPtr);
895 template<
class V,
class M>
897 BaseVectorSequence<V,M>::computeMeanVars(
898 const SsOptionsValues& statisticalOptions,
899 std::ofstream* passedOfs,
906 struct timeval timevalTmp;
907 iRC = gettimeofday(&timevalTmp, NULL);
908 double tmpRunTime = 0.;
910 if (m_env.subDisplayFile()) {
911 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
912 <<
"\nComputing mean, sample variance and population variance"
916 V subChainMean(m_vectorSpace.zeroVector());
917 this->subMeanExtra(0,
918 this->subSequenceSize(),
921 V subChainMedian(m_vectorSpace.zeroVector());
922 this->subMedianExtra(0,
923 this->subSequenceSize(),
926 V subChainSampleVariance(m_vectorSpace.zeroVector());
927 this->subSampleVarianceExtra(0,
928 this->subSequenceSize(),
930 subChainSampleVariance);
932 if ((m_env.displayVerbosity() >= 5) && (m_env.subDisplayFile())) {
933 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeMeanVars()"
934 <<
": subChainMean.sizeLocal() = " << subChainMean.sizeLocal()
935 <<
", subChainMean = " << subChainMean
936 <<
", subChainMedian = " << subChainMedian
937 <<
", subChainSampleVariance.sizeLocal() = " << subChainSampleVariance.sizeLocal()
938 <<
", subChainSampleVariance = " << subChainSampleVariance
942 V estimatedVarianceOfSampleMean(subChainSampleVariance);
943 estimatedVarianceOfSampleMean /= (double) this->subSequenceSize();
944 bool savedVectorPrintState = estimatedVarianceOfSampleMean.getPrintHorizontally();
945 estimatedVarianceOfSampleMean.setPrintHorizontally(
false);
946 if (m_env.subDisplayFile()) {
947 *m_env.subDisplayFile() <<
"\nEstimated variance of sample mean for the whole chain '" << m_name <<
"'"
948 <<
", under independence assumption:"
950 << estimatedVarianceOfSampleMean
953 estimatedVarianceOfSampleMean.setPrintHorizontally(savedVectorPrintState);
955 V estimatedStdOfSampleMean(estimatedVarianceOfSampleMean);
956 estimatedStdOfSampleMean.cwSqrt();
957 savedVectorPrintState = estimatedStdOfSampleMean.getPrintHorizontally();
958 estimatedStdOfSampleMean.setPrintHorizontally(
false);
959 if (m_env.subDisplayFile()) {
960 *m_env.subDisplayFile() <<
"\nEstimated standard deviation of sample mean for the whole chain '" << m_name <<
"'"
961 <<
", under independence assumption:"
963 << estimatedStdOfSampleMean
966 estimatedStdOfSampleMean.setPrintHorizontally(savedVectorPrintState);
968 V subChainPopulationVariance(m_vectorSpace.zeroVector());
969 this->subPopulationVariance(0,
970 this->subSequenceSize(),
972 subChainPopulationVariance);
975 if (m_env.subDisplayFile()) {
976 *m_env.subDisplayFile() <<
"Sub Mean, median, and variances took " << tmpRunTime
981 if (m_env.subDisplayFile()) {
982 *m_env.subDisplayFile() <<
"\nSub mean, median, sample std, population std"
985 sprintf(line,
"%s%4s%s%6s%s%9s%s%9s%s",
995 *m_env.subDisplayFile() << line;
997 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
998 sprintf(line,
"\n%8.8s%2s%11.4e%2s%11.4e%2s%11.4e%2s%11.4e",
999 m_vectorSpace.localComponentName(i).c_str(),
1005 std::sqrt(subChainSampleVariance[i]),
1007 std::sqrt(subChainPopulationVariance[i]));
1008 *m_env.subDisplayFile() << line;
1010 *m_env.subDisplayFile() << std::endl;
1013 if (subMeanPtr ) *subMeanPtr = subChainMean;
1014 if (subMedianPtr ) *subMedianPtr = subChainMedian;
1015 if (subSampleVarPtr) *subSampleVarPtr = subChainSampleVariance;
1016 if (subPopulVarPtr ) *subPopulVarPtr = subChainPopulationVariance;
1018 if (m_env.numSubEnvironments() > 1) {
1020 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1021 V unifiedChainMean(m_vectorSpace.zeroVector());
1022 this->unifiedMeanExtra(0,
1023 this->subSequenceSize(),
1026 V unifiedChainMedian(m_vectorSpace.zeroVector());
1027 this->unifiedMedianExtra(0,
1028 this->subSequenceSize(),
1029 unifiedChainMedian);
1031 V unifiedChainSampleVariance(m_vectorSpace.zeroVector());
1032 this->unifiedSampleVarianceExtra(0,
1033 this->subSequenceSize(),
1035 unifiedChainSampleVariance);
1037 V unifiedChainPopulationVariance(m_vectorSpace.zeroVector());
1038 this->unifiedPopulationVariance(0,
1039 this->subSequenceSize(),
1041 unifiedChainPopulationVariance);
1043 if (m_env.inter0Rank() == 0) {
1044 if (m_env.subDisplayFile()) {
1045 *m_env.subDisplayFile() <<
"\nUnif mean, median, sample std, population std"
1048 sprintf(line,
"%s%4s%s%6s%s%9s%s%9s%s",
1058 *m_env.subDisplayFile() << line;
1060 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1061 sprintf(line,
"\n%8.8s%2s%11.4e%2s%11.4e%2s%11.4e%2s%11.4e",
1062 m_vectorSpace.localComponentName(i).c_str(),
1064 unifiedChainMean[i],
1066 unifiedChainMedian[i],
1068 std::sqrt(unifiedChainSampleVariance[i]),
1070 std::sqrt(unifiedChainPopulationVariance[i]));
1071 *m_env.subDisplayFile() << line;
1073 *m_env.subDisplayFile() << std::endl;
1078 queso_error_msg(
"unified min-max writing, parallel vectors not supported yet");
1082 if (m_env.subDisplayFile()) {
1083 *m_env.subDisplayFile() <<
"\nEnded computing mean, sample variance and population variance"
1084 <<
"\n-----------------------------------------------------"
1091 template<
class V,
class M>
1093 BaseVectorSequence<V,M>::computeAutoCorrViaDef(
1094 const SequenceStatisticalOptions& statisticalOptions,
1095 const std::vector<unsigned int>& initialPosForStatistics,
1096 const std::vector<unsigned int>& lagsForCorrs,
1097 std::ofstream* passedOfs)
1100 computeAutoCorrViaDef(statisticalOptions.m_ov, initialPosForStatistics,
1101 lagsForCorrs, passedOfs);
1104 template<
class V,
class M>
1106 BaseVectorSequence<V,M>::computeAutoCorrViaDef(
1107 const SsOptionsValues& statisticalOptions,
1108 const std::vector<unsigned int>& initialPosForStatistics,
1109 const std::vector<unsigned int>& lagsForCorrs,
1110 std::ofstream* passedOfs)
1113 struct timeval timevalTmp;
1114 iRC = gettimeofday(&timevalTmp, NULL);
1115 double tmpRunTime = 0.;
1117 if (m_env.subDisplayFile()) {
1118 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1119 <<
"\nComputing autocorrelation coefficients (via def)"
1123 if (statisticalOptions.autoCorrDisplay() && (m_env.subDisplayFile())) {
1124 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeAutoCorrViaDef(): lags for autocorrelation (via def) = ";
1125 for (
unsigned int i = 0; i < lagsForCorrs.size(); ++i) {
1126 *m_env.subDisplayFile() <<
" " << lagsForCorrs[i];
1128 *m_env.subDisplayFile() << std::endl;
1131 TwoDArray<V> _2dArrayOfAutoCorrs(initialPosForStatistics.size(),lagsForCorrs.size());
1132 for (
unsigned int i = 0; i < _2dArrayOfAutoCorrs.numRows(); ++i) {
1133 for (
unsigned int j = 0; j < _2dArrayOfAutoCorrs.numCols(); ++j) {
1134 _2dArrayOfAutoCorrs.setLocation(i,j,
new V(m_vectorSpace.zeroVector()) );
1138 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1139 unsigned int initialPos = initialPosForStatistics[initialPosId];
1140 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1141 unsigned int lag = lagsForCorrs[lagId];
1142 this->autoCorrViaDef(initialPos,
1143 this->subSequenceSize()-initialPos,
1145 _2dArrayOfAutoCorrs(initialPosId,lagId));
1153 if ((statisticalOptions.autoCorrDisplay()) && (m_env.subDisplayFile())) {
1154 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1155 *m_env.subDisplayFile() <<
"\nComputed autocorrelation coefficients (via def), for subchain beginning at position " << initialPosForStatistics[initialPosId]
1156 <<
" (each column corresponds to a different lag)"
1161 *m_env.subDisplayFile() << line;
1162 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1163 sprintf(line,
"%10s%3d",
1165 lagsForCorrs[lagId]);
1166 *m_env.subDisplayFile() << line;
1169 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1170 sprintf(line,
"\n%9.9s",
1171 m_vectorSpace.localComponentName(i).c_str() );
1172 *m_env.subDisplayFile() << line;
1173 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1174 sprintf(line,
"%2s%11.4e",
1176 _2dArrayOfAutoCorrs(initialPosId,lagId)[i]);
1177 *m_env.subDisplayFile() << line;
1180 *m_env.subDisplayFile() << std::endl;
1185 if (m_env.subDisplayFile()) {
1186 *m_env.subDisplayFile() <<
"Chain autocorrelation (via def) took " << tmpRunTime
1192 if (statisticalOptions.autoCorrWrite() && passedOfs) {
1193 std::ofstream& ofsvar = *passedOfs;
1194 ofsvar << m_name <<
"_corrViaDefLags_sub" << m_env.subIdString() <<
" = zeros(" << 1
1195 <<
"," << lagsForCorrs.size()
1198 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1199 ofsvar << m_name <<
"_corrViaDefLags_sub" << m_env.subIdString() <<
"(" << 1
1201 <<
") = " << lagsForCorrs[lagId]
1206 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1207 ofsvar << m_name <<
"_corrViaDefInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << this->vectorSizeLocal()
1208 <<
"," << lagsForCorrs.size()
1211 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1212 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1213 ofsvar << m_name <<
"_corrViaDefInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << i+1
1215 <<
") = " << _2dArrayOfAutoCorrs(initialPosId,lagId)[i]
1226 template<
class V,
class M>
1228 BaseVectorSequence<V,M>::computeAutoCorrViaFFT(
1229 const SequenceStatisticalOptions& statisticalOptions,
1230 const std::vector<unsigned int>& initialPosForStatistics,
1231 const std::vector<unsigned int>& lagsForCorrs,
1232 std::ofstream* passedOfs)
1235 computeAutoCorrViaFFT(statisticalOptions.m_ov, initialPosForStatistics,
1236 lagsForCorrs, passedOfs);
1239 template<
class V,
class M>
1241 BaseVectorSequence<V,M>::computeAutoCorrViaFFT(
1242 const SsOptionsValues& statisticalOptions,
1243 const std::vector<unsigned int>& initialPosForStatistics,
1244 const std::vector<unsigned int>& lagsForCorrs,
1245 std::ofstream* passedOfs)
1248 struct timeval timevalTmp;
1249 iRC = gettimeofday(&timevalTmp, NULL);
1250 double tmpRunTime = 0.;
1252 if (m_env.subDisplayFile()) {
1253 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1254 <<
"\nComputing autocorrelation coefficients (via fft)"
1258 if (statisticalOptions.autoCorrDisplay() && (m_env.subDisplayFile())) {
1259 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeAutoCorrViaFFT(): lags for autocorrelation (via fft) = ";
1260 for (
unsigned int i = 0; i < lagsForCorrs.size(); ++i) {
1261 *m_env.subDisplayFile() <<
" " << lagsForCorrs[i];
1263 *m_env.subDisplayFile() << std::endl;
1266 TwoDArray<V> _2dArrayOfAutoCorrs(initialPosForStatistics.size(),lagsForCorrs.size());
1267 for (
unsigned int i = 0; i < _2dArrayOfAutoCorrs.numRows(); ++i) {
1268 for (
unsigned int j = 0; j < _2dArrayOfAutoCorrs.numCols(); ++j) {
1269 _2dArrayOfAutoCorrs.setLocation(i,j,
new V(m_vectorSpace.zeroVector()) );
1272 std::vector<V*> corrVecs(lagsForCorrs.size(),NULL);
1273 std::vector<V*> corrSumVecs(initialPosForStatistics.size(),NULL);
1274 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1275 corrSumVecs[initialPosId] =
new V(m_vectorSpace.zeroVector()) ;
1276 unsigned int initialPos = initialPosForStatistics[initialPosId];
1277 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1278 corrVecs[lagId] =
new V(m_vectorSpace.zeroVector()) ;
1280 if (m_env.subDisplayFile()) {
1281 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeAutoCorrViaFFT()"
1282 <<
": about to call chain.autoCorrViaFft()"
1283 <<
" with initialPos = " << initialPos
1284 <<
", numPos = " << this->subSequenceSize()-initialPos
1285 <<
", lagsForCorrs.size() = " << lagsForCorrs.size()
1286 <<
", corrVecs.size() = " << corrVecs.size()
1289 this->autoCorrViaFft(initialPos,
1290 this->subSequenceSize()-initialPos,
1293 this->autoCorrViaFft(initialPos,
1294 this->subSequenceSize()-initialPos,
1295 (
unsigned int) (1.0 * (
double) (this->subSequenceSize()-initialPos)),
1296 *corrSumVecs[initialPosId]);
1297 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1298 _2dArrayOfAutoCorrs(initialPosId,lagId) = *(corrVecs[lagId]);
1301 for (
unsigned int j = 0; j < corrVecs.size(); ++j) {
1302 if (corrVecs[j] != NULL)
delete corrVecs[j];
1305 if (statisticalOptions.autoCorrDisplay()) {
1306 V subChainMean (m_vectorSpace.zeroVector());
1307 V subChainSampleVariance (m_vectorSpace.zeroVector());
1308 V estimatedVarianceOfSampleMean(m_vectorSpace.zeroVector());
1309 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1310 unsigned int initialPos = initialPosForStatistics[initialPosId];
1312 this->subMeanExtra(initialPos,
1313 this->subSequenceSize()-initialPos,
1316 this->subSampleVarianceExtra(initialPos,
1317 this->subSequenceSize()-initialPos,
1319 subChainSampleVariance);
1321 if (m_env.subDisplayFile()) {
1322 *m_env.subDisplayFile() <<
"\nEstimated variance of sample mean, through autocorrelation (via fft), for subchain beginning at position " << initialPosForStatistics[initialPosId]
1325 estimatedVarianceOfSampleMean.cwSet(-1.);
1326 estimatedVarianceOfSampleMean += 2.* (*corrSumVecs[initialPosId]);
1327 estimatedVarianceOfSampleMean *= subChainSampleVariance;
1328 estimatedVarianceOfSampleMean /= (double) (this->subSequenceSize() - initialPos);
1329 bool savedVectorPrintState = estimatedVarianceOfSampleMean.getPrintHorizontally();
1330 estimatedVarianceOfSampleMean.setPrintHorizontally(
false);
1331 if (m_env.subDisplayFile()) {
1332 *m_env.subDisplayFile() << estimatedVarianceOfSampleMean
1335 estimatedVarianceOfSampleMean.setPrintHorizontally(savedVectorPrintState);
1337 if (m_env.subDisplayFile()) {
1338 *m_env.subDisplayFile() <<
"\nComputed autocorrelation coefficients (via fft), for subchain beginning at position " << initialPosForStatistics[initialPosId]
1339 <<
" (each column corresponds to a different lag)"
1345 *m_env.subDisplayFile() << line;
1346 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1347 sprintf(line,
"%10s%3d",
1349 lagsForCorrs[lagId]);
1350 *m_env.subDisplayFile() << line;
1353 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1354 sprintf(line,
"\n%9.9s",
1355 m_vectorSpace.localComponentName(i).c_str() );
1356 *m_env.subDisplayFile() << line;
1357 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1358 sprintf(line,
"%2s%11.4e",
1360 _2dArrayOfAutoCorrs(initialPosId,lagId)[i]);
1361 *m_env.subDisplayFile() << line;
1364 *m_env.subDisplayFile() << std::endl;
1370 if (m_env.subDisplayFile()) {
1371 *m_env.subDisplayFile() <<
"Chain autocorrelation (via fft) took " << tmpRunTime
1377 if (statisticalOptions.autoCorrWrite() && passedOfs) {
1378 std::ofstream& ofsvar = *passedOfs;
1379 ofsvar << m_name <<
"_corrViaFftLags_sub" << m_env.subIdString() <<
" = zeros(" << 1
1380 <<
"," << lagsForCorrs.size()
1383 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1384 ofsvar << m_name <<
"_corrViaFftLags_sub" << m_env.subIdString() <<
"(" << 1
1386 <<
") = " << lagsForCorrs[lagId]
1391 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1392 ofsvar << m_name <<
"_corrViaFftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << this->vectorSizeLocal()
1393 <<
"," << lagsForCorrs.size()
1396 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1397 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1398 ofsvar << m_name <<
"_corrViaFftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << i+1
1400 <<
") = " << _2dArrayOfAutoCorrs(initialPosId,lagId)[i]
1411 template<
class V,
class M>
1413 BaseVectorSequence<V,M>::computeHistCdfstaccKde(
1414 const SequenceStatisticalOptions& statisticalOptions,
1415 std::ofstream* passedOfs)
1418 computeHistCdfstaccKde(statisticalOptions.m_ov, passedOfs);
1421 template<
class V,
class M>
1423 BaseVectorSequence<V,M>::computeHistCdfstaccKde(
1424 const SsOptionsValues& statisticalOptions,
1425 std::ofstream* passedOfs)
1427 if (m_env.subDisplayFile()) {
1428 *m_env.subDisplayFile() <<
"\n"
1429 <<
"\n-----------------------------------------------------"
1430 <<
"\n Computing histogram and/or cdf stacc and/or Kde for chain '" << m_name <<
"' ..."
1431 <<
"\n-----------------------------------------------------"
1437 struct timeval timevalTmp;
1442 double tmpRunTime = 0.;
1443 iRC = gettimeofday(&timevalTmp, NULL);
1444 if (m_env.subDisplayFile()) {
1445 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1446 <<
"\nComputing min and max for histograms and Kde"
1450 V statsMinPositions(m_vectorSpace.zeroVector());
1451 V statsMaxPositions(m_vectorSpace.zeroVector());
1452 this->subMinMaxExtra(0,
1453 this->subSequenceSize(),
1457 if (m_env.subDisplayFile()) {
1458 *m_env.subDisplayFile() <<
"\nComputed min values and max values for chain '" << m_name <<
"'"
1464 *m_env.subDisplayFile() << line;
1466 sprintf(line,
"%9s%s%9s%s",
1471 *m_env.subDisplayFile() << line;
1473 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1474 sprintf(line,
"\n%8.8s",
1475 m_vectorSpace.localComponentName(i).c_str() );
1476 *m_env.subDisplayFile() << line;
1478 sprintf(line,
"%2s%11.4e%2s%11.4e",
1480 statsMinPositions[i],
1482 statsMaxPositions[i]);
1483 *m_env.subDisplayFile() << line;
1485 *m_env.subDisplayFile() << std::endl;
1488 V unifiedStatsMinPositions(statsMinPositions);
1489 V unifiedStatsMaxPositions(statsMaxPositions);
1490 if (m_env.numSubEnvironments() > 1) {
1492 this->unifiedMinMaxExtra(0,
1493 this->subSequenceSize(),
1494 unifiedStatsMinPositions,
1495 unifiedStatsMaxPositions);
1498 if (m_env.subDisplayFile()) {
1499 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1500 if (m_env.inter0Rank() == 0) {
1501 *m_env.subDisplayFile() <<
"\nComputed unified min values and max values for chain '" << m_name <<
"'"
1507 *m_env.subDisplayFile() << line;
1509 sprintf(line,
"%9s%s%9s%s",
1514 *m_env.subDisplayFile() << line;
1516 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1517 sprintf(line,
"\n%8.8s",
1518 m_vectorSpace.localComponentName(i).c_str() );
1519 *m_env.subDisplayFile() << line;
1521 sprintf(line,
"%2s%11.4e%2s%11.4e",
1523 unifiedStatsMinPositions[i],
1525 unifiedStatsMaxPositions[i]);
1526 *m_env.subDisplayFile() << line;
1528 *m_env.subDisplayFile() << std::endl;
1532 queso_error_msg(
"unified min-max writing, parallel vectors not supported yet");
1539 if (m_env.subDisplayFile()) {
1540 *m_env.subDisplayFile() <<
"Chain min and max took " << tmpRunTime
1548 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
1549 if ((statisticalOptions.histCompute() ) &&
1550 (statisticalOptions.histNumInternalBins() > 0)) {
1552 iRC = gettimeofday(&timevalTmp, NULL);
1553 if (m_env.subDisplayFile()) {
1554 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1555 <<
"\nComputing histograms"
1559 std::string subCoreName_HistCenters((std::string)(
"_HistCenters_sub")+m_env.subIdString());
1560 std::string uniCoreName_HistCenters((std::string)(
"_unifHistCenters_sub")+m_env.subIdString());
1561 if (m_env.numSubEnvironments() == 1) subCoreName_HistCenters = uniCoreName_HistCenters;
1563 std::string subCoreName_HistQuantts((std::string)(
"_HistQuantts_sub")+m_env.subIdString());
1564 std::string uniCoreName_HistQuantts((std::string)(
"_unifHistQuantts_sub")+m_env.subIdString());
1565 if (m_env.numSubEnvironments() == 1) subCoreName_HistQuantts = uniCoreName_HistQuantts;
1567 for (
unsigned int i = 0; i < statsMaxPositions.sizeLocal(); ++i) {
1568 statsMaxPositions[i] *= (1. + 1.e-15);
1572 std::vector<V*> histCentersForAllBins(statisticalOptions.histNumInternalBins()+2,NULL);
1573 std::vector<V*> histQuanttsForAllBins(statisticalOptions.histNumInternalBins()+2,NULL);
1574 this->subHistogram(0,
1577 histCentersForAllBins,
1578 histQuanttsForAllBins);
1582 std::ofstream& ofsvar = *passedOfs;
1583 ofsvar << m_name << subCoreName_HistCenters <<
" = zeros(" << this->vectorSizeLocal()
1584 <<
"," << histCentersForAllBins.size()
1587 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1588 for (
unsigned int j = 0; j < histCentersForAllBins.size(); ++j) {
1589 ofsvar << m_name << subCoreName_HistCenters <<
"(" << i+1
1591 <<
") = " << (*(histCentersForAllBins[j]))[i]
1597 ofsvar << m_name << subCoreName_HistQuantts <<
" = zeros(" << this->vectorSizeLocal()
1598 <<
"," << histQuanttsForAllBins.size()
1601 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1602 for (
unsigned int j = 0; j < histQuanttsForAllBins.size(); ++j) {
1603 ofsvar << m_name << subCoreName_HistQuantts <<
"(" << i+1
1605 <<
") = " << (*(histQuanttsForAllBins[j]))[i]
1612 for (
unsigned int i = 0; i < histQuanttsForAllBins.size(); ++i) {
1613 if (histQuanttsForAllBins[i] != NULL)
delete histQuanttsForAllBins[i];
1615 for (
unsigned int i = 0; i < histCentersForAllBins.size(); ++i) {
1616 if (histCentersForAllBins[i] != NULL)
delete histCentersForAllBins[i];
1619 std::vector<V*> unifiedHistCentersForAllBins(statisticalOptions.histNumInternalBins()+2,NULL);
1620 std::vector<V*> unifiedHistQuanttsForAllBins(statisticalOptions.histNumInternalBins()+2,NULL);
1621 if (m_env.numSubEnvironments() > 1) {
1623 this->unifiedHistogram(0,
1624 unifiedStatsMinPositions,
1625 unifiedStatsMaxPositions,
1626 unifiedHistCentersForAllBins,
1627 unifiedHistQuanttsForAllBins);
1631 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1632 if (m_env.inter0Rank() == 0) {
1633 std::ofstream& ofsvar = *passedOfs;
1634 ofsvar << m_name << uniCoreName_HistCenters <<
" = zeros(" << this->vectorSizeLocal()
1635 <<
"," << unifiedHistCentersForAllBins.size()
1638 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1639 for (
unsigned int j = 0; j < unifiedHistCentersForAllBins.size(); ++j) {
1640 ofsvar << m_name << uniCoreName_HistCenters <<
"(" << i+1
1642 <<
") = " << (*(unifiedHistCentersForAllBins[j]))[i]
1648 ofsvar << m_name << uniCoreName_HistQuantts <<
" = zeros(" << this->vectorSizeLocal()
1649 <<
"," << unifiedHistQuanttsForAllBins.size()
1652 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1653 for (
unsigned int j = 0; j < unifiedHistQuanttsForAllBins.size(); ++j) {
1654 ofsvar << m_name << uniCoreName_HistQuantts <<
"(" << i+1
1656 <<
") = " << (*(unifiedHistQuanttsForAllBins[j]))[i]
1664 queso_error_msg(
"unified histogram writing, parallel vectors not supported yet");
1668 for (
unsigned int i = 0; i < unifiedHistQuanttsForAllBins.size(); ++i) {
1669 if (unifiedHistQuanttsForAllBins[i] != NULL)
delete unifiedHistQuanttsForAllBins[i];
1671 for (
unsigned int i = 0; i < unifiedHistCentersForAllBins.size(); ++i) {
1672 if (unifiedHistCentersForAllBins[i] != NULL)
delete unifiedHistCentersForAllBins[i];
1678 if (m_env.subDisplayFile()) {
1679 *m_env.subDisplayFile() <<
"Chain histograms took " << tmpRunTime
1685 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
1689 if ((statisticalOptions.cdfStaccCompute() ) &&
1690 (statisticalOptions.cdfStaccNumEvalPositions() > 0)) {
1692 iRC = gettimeofday(&timevalTmp, NULL);
1693 if (m_env.subDisplayFile()) {
1694 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1695 <<
"\nComputing cdf statistical accuracy"
1699 std::vector<V*> cdfStaccEvalPositions(statisticalOptions.cdfStaccNumEvalPositions(),NULL);
1702 cdfStaccEvalPositions);
1704 std::vector<V*> cdfStaccValues(statisticalOptions.cdfStaccNumEvalPositions(),NULL);
1705 this->subCdfStacc(0,
1706 cdfStaccEvalPositions,
1711 if (m_env.subDisplayFile()) {
1712 *m_env.subDisplayFile() <<
"Chain cdf statistical accuracy took " << tmpRunTime
1721 if ((statisticalOptions.kdeCompute() ) &&
1722 (statisticalOptions.kdeNumEvalPositions() > 0)) {
1724 iRC = gettimeofday(&timevalTmp, NULL);
1725 if (m_env.subDisplayFile()) {
1726 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1727 <<
"\nComputing Kde"
1731 std::string subCoreName_GaussianKdePositions((std::string)(
"_GkdePosits_sub")+m_env.subIdString());
1732 std::string uniCoreName_GaussianKdePositions((std::string)(
"_unifGkdePosits_sub")+m_env.subIdString());
1733 if (m_env.numSubEnvironments() == 1) subCoreName_GaussianKdePositions = uniCoreName_GaussianKdePositions;
1735 std::string subCoreName_GaussianKdeScaleVec ((std::string)(
"_GkdeScalev_sub")+m_env.subIdString());
1736 std::string uniCoreName_GaussianKdeScaleVec ((std::string)(
"_unifGkdeScalev_sub")+m_env.subIdString());
1737 if (m_env.numSubEnvironments() == 1) subCoreName_GaussianKdeScaleVec = uniCoreName_GaussianKdeScaleVec;
1739 std::string subCoreName_GaussianKdeValues ((std::string)(
"_GkdeValues_sub")+m_env.subIdString());
1740 std::string uniCoreName_GaussianKdeValues ((std::string)(
"_unifGkdeValues_sub")+m_env.subIdString());
1741 if (m_env.numSubEnvironments() == 1) subCoreName_GaussianKdeValues = uniCoreName_GaussianKdeValues;
1743 V iqrVec(m_vectorSpace.zeroVector());
1744 this->subInterQuantileRange(0,
1747 V gaussianKdeScaleVec(m_vectorSpace.zeroVector());
1748 this->subScalesForKde(0,
1751 gaussianKdeScaleVec);
1753 std::vector<V*> kdeEvalPositions(statisticalOptions.kdeNumEvalPositions(),NULL);
1758 std::vector<V*> gaussianKdeDensities(statisticalOptions.kdeNumEvalPositions(),NULL);
1759 this->subGaussian1dKde(0,
1760 gaussianKdeScaleVec,
1762 gaussianKdeDensities);
1765 if (m_env.subDisplayFile()) {
1766 *m_env.subDisplayFile() <<
"\nComputed inter quantile ranges for chain '" << m_name <<
"'"
1772 *m_env.subDisplayFile() << line;
1774 sprintf(line,
"%9s%s",
1777 *m_env.subDisplayFile() << line;
1779 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1780 sprintf(line,
"\n%8.8s",
1781 m_vectorSpace.localComponentName(i).c_str() );
1782 *m_env.subDisplayFile() << line;
1784 sprintf(line,
"%2s%11.4e",
1787 *m_env.subDisplayFile() << line;
1789 *m_env.subDisplayFile() << std::endl;
1793 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() > 10)) {
1794 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeHistCdfstaccKde()"
1795 <<
", for chain '" << m_name <<
"'"
1796 <<
", about to write sub kde to ofstream with pointer = " << passedOfs
1800 std::ofstream& ofsvar = *passedOfs;
1801 ofsvar << m_name << subCoreName_GaussianKdePositions <<
" = zeros(" << this->vectorSizeLocal()
1802 <<
"," << kdeEvalPositions.size()
1805 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1806 for (
unsigned int j = 0; j < kdeEvalPositions.size(); ++j) {
1807 ofsvar << m_name << subCoreName_GaussianKdePositions <<
"(" << i+1
1809 <<
") = " << (*(kdeEvalPositions[j]))[i]
1815 ofsvar << m_name << subCoreName_GaussianKdeScaleVec <<
" = zeros(" << this->vectorSizeLocal()
1818 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1819 ofsvar << m_name << subCoreName_GaussianKdeScaleVec <<
"(" << i+1
1820 <<
") = " << gaussianKdeScaleVec[i]
1825 ofsvar << m_name << subCoreName_GaussianKdeValues <<
" = zeros(" << this->vectorSizeLocal()
1826 <<
"," << gaussianKdeDensities.size()
1829 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1830 for (
unsigned int j = 0; j < gaussianKdeDensities.size(); ++j) {
1831 ofsvar << m_name << subCoreName_GaussianKdeValues <<
"(" << i+1
1833 <<
") = " << (*(gaussianKdeDensities[j]))[i]
1840 for (
unsigned int i = 0; i < gaussianKdeDensities.size(); ++i) {
1841 if (gaussianKdeDensities[i] != NULL)
delete gaussianKdeDensities[i];
1843 for (
unsigned int i = 0; i < kdeEvalPositions.size(); ++i) {
1844 if (kdeEvalPositions[i] != NULL)
delete kdeEvalPositions[i];
1847 if ((
int) m_env.numSubEnvironments() > 1) {
1849 V unifiedIqrVec(m_vectorSpace.zeroVector());
1850 this->unifiedInterQuantileRange(0,
1854 V unifiedGaussianKdeScaleVec(m_vectorSpace.zeroVector());
1855 this->unifiedScalesForKde(0,
1858 unifiedGaussianKdeScaleVec);
1861 std::vector<V*> unifiedKdeEvalPositions(statisticalOptions.kdeNumEvalPositions(),NULL);
1863 unifiedStatsMaxPositions,
1864 unifiedKdeEvalPositions);
1866 std::vector<V*> unifiedGaussianKdeDensities(statisticalOptions.kdeNumEvalPositions(),NULL);
1867 this->unifiedGaussian1dKde(0,
1868 unifiedGaussianKdeScaleVec,
1869 unifiedKdeEvalPositions,
1870 unifiedGaussianKdeDensities);
1874 if (m_env.subDisplayFile()) {
1875 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1876 if (m_env.inter0Rank() == 0) {
1877 *m_env.subDisplayFile() <<
"\nComputed unified inter quantile ranges for chain '" << m_name <<
"'"
1883 *m_env.subDisplayFile() << line;
1885 sprintf(line,
"%9s%s",
1888 *m_env.subDisplayFile() << line;
1890 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1891 sprintf(line,
"\n%8.8s",
1892 m_vectorSpace.localComponentName(i).c_str() );
1893 *m_env.subDisplayFile() << line;
1895 sprintf(line,
"%2s%11.4e",
1898 *m_env.subDisplayFile() << line;
1900 *m_env.subDisplayFile() << std::endl;
1904 queso_error_msg(
"unified iqr writing, parallel vectors not supported yet");
1909 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() > 10)) {
1910 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeHistCdfstaccKde()"
1911 <<
", for chain '" << m_name <<
"'"
1912 <<
", about to write unified kde to ofstream with pointer = " << passedOfs
1916 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1917 if (m_env.inter0Rank() == 0) {
1918 std::ofstream& ofsvar = *passedOfs;
1919 ofsvar << m_name << uniCoreName_GaussianKdePositions <<
" = zeros(" << this->vectorSizeLocal()
1920 <<
"," << unifiedKdeEvalPositions.size()
1923 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() > 10)) {
1924 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeHistCdfstaccKde()"
1925 <<
", for chain '" << m_name <<
"'"
1926 <<
": just wrote '... = zeros(.,.);' line to output file, which has pointer = " << passedOfs
1929 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1930 for (
unsigned int j = 0; j < unifiedKdeEvalPositions.size(); ++j) {
1931 ofsvar << m_name << uniCoreName_GaussianKdePositions <<
"(" << i+1
1933 <<
") = " << (*(unifiedKdeEvalPositions[j]))[i]
1939 ofsvar << m_name << uniCoreName_GaussianKdeScaleVec <<
" = zeros(" << this->vectorSizeLocal()
1942 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1943 ofsvar << m_name << uniCoreName_GaussianKdeScaleVec <<
"(" << i+1
1944 <<
") = " << unifiedGaussianKdeScaleVec[i]
1949 ofsvar << m_name << uniCoreName_GaussianKdeValues <<
" = zeros(" << this->vectorSizeLocal()
1950 <<
"," << unifiedGaussianKdeDensities.size()
1953 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1954 for (
unsigned int j = 0; j < unifiedGaussianKdeDensities.size(); ++j) {
1955 ofsvar << m_name << uniCoreName_GaussianKdeValues <<
"(" << i+1
1957 <<
") = " << (*(unifiedGaussianKdeDensities[j]))[i]
1965 queso_error_msg(
"unified Kde writing, parallel vectors not supported yet");
1969 for (
unsigned int i = 0; i < unifiedGaussianKdeDensities.size(); ++i) {
1970 if (unifiedGaussianKdeDensities[i] != NULL)
delete unifiedGaussianKdeDensities[i];
1972 for (
unsigned int i = 0; i < unifiedKdeEvalPositions.size(); ++i) {
1973 if (unifiedKdeEvalPositions[i] != NULL)
delete unifiedKdeEvalPositions[i];
1979 if (m_env.subDisplayFile()) {
1980 *m_env.subDisplayFile() <<
"Chain Kde took " << tmpRunTime
1986 if (m_env.subDisplayFile()) {
1987 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1988 <<
"\n Finished computing histogram and/or cdf stacc and/or Kde for chain '" << m_name <<
"'"
1989 <<
"\n-----------------------------------------------------"
1998 template<
class V,
class M>
2000 BaseVectorSequence<V,M>::computeCovCorrMatrices(
2001 const SequenceStatisticalOptions& statisticalOptions,
2002 std::ofstream* passedOfs)
2005 computeCovCorrMatrices(statisticalOptions.m_ov, passedOfs);
2008 template<
class V,
class M>
2010 BaseVectorSequence<V,M>::computeCovCorrMatrices(
2011 const SsOptionsValues& statisticalOptions,
2012 std::ofstream* passedOfs)
2014 if (m_env.subDisplayFile()) {
2015 *m_env.subDisplayFile() <<
"\n"
2016 <<
"\n-----------------------------------------------------"
2017 <<
"\n Computing covariance and correlation matrices for chain '" << m_name <<
"' ..."
2018 <<
"\n-----------------------------------------------------"
2025 M* covarianceMatrix =
new M(m_env,
2026 m_vectorSpace.map(),
2027 m_vectorSpace.dimGlobal());
2028 M* correlationMatrix =
new M(m_env,
2029 m_vectorSpace.map(),
2030 m_vectorSpace.dimGlobal());
2034 this->subSequenceSize(),
2036 *correlationMatrix);
2038 if (m_env.subDisplayFile()) {
2039 if (m_vectorSpace.numOfProcsForStorage() == 1) {
2041 if (m_env.inter0Rank() == 0) {
2042 *m_env.subDisplayFile() <<
"\nBaseVectorSequence<V,M>::computeCovCorrMatrices"
2043 <<
", chain " << m_name
2044 <<
": contents of covariance matrix are\n" << *covarianceMatrix
2047 *m_env.subDisplayFile() <<
"\nBaseVectorSequence<V,M>::computeCovCorrMatrices"
2048 <<
", chain " << m_name
2049 <<
": contents of correlation matrix are\n" << *correlationMatrix
2058 delete correlationMatrix;
2059 delete covarianceMatrix;
2061 if (m_env.subDisplayFile()) {
2062 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2063 <<
"\n Finished computing covariance and correlation matrices for chain '" << m_name <<
"'"
2064 <<
"\n-----------------------------------------------------"
2071 #endif // #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2077 #ifdef UQ_CODE_HAS_MONITORS
2078 template<
class V,
class M>
2080 BaseVectorSequence<V,M>::computeMeanEvolution(
2081 const SequenceStatisticalOptions& statisticalOptions,
2082 std::ofstream* passedOfs)
2085 computeMeanEvolution(statisticalOptions.m_ov, passedOfs);
2088 template<
class V,
class M>
2090 BaseVectorSequence<V,M>::computeMeanEvolution(
2091 const SsOptionsValues& statisticalOptions,
2092 std::ofstream* passedOfs)
2095 struct timeval timevalTmp;
2096 iRC = gettimeofday(&timevalTmp, NULL);
2097 double tmpRunTime = 0.;
2099 if (m_env.subDisplayFile()) {
2100 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2101 <<
"\nComputing mean evolution"
2105 unsigned int monitorPeriod = statisticalOptions.meanMonitorPeriod();
2106 unsigned int iMin = 0;
2107 unsigned int iMax = (
unsigned int) ( ((
double) this->subSequenceSize())/((double) monitorPeriod) );
2109 if (monitorPeriod == 1) {
2113 if (m_env.subDisplayFile()) {
2114 *m_env.subDisplayFile() <<
" Sub sequence size = " << this->subSequenceSize()
2115 <<
"\n Monitor period = " << monitorPeriod
2116 <<
"\n Number of monitoring positions = " << iMax
2120 this->subMeanMonitorAlloc(iMax);
2121 if (m_env.numSubEnvironments() > 1) {
2122 this->subMeanInter0MonitorAlloc(iMax);
2123 this->unifiedMeanMonitorAlloc(iMax);
2126 for (
unsigned int i = iMin; i < iMax; ++i) {
2127 unsigned int currentMonitoredFinalPosition = (i+1)*monitorPeriod - 1;
2128 V subMeanVec (m_vectorSpace.zeroVector());
2129 V subMeanCltStd(m_vectorSpace.zeroVector());
2130 this->subMeanMonitorRun(currentMonitoredFinalPosition,
2133 this->subMeanMonitorStore(i,
2134 currentMonitoredFinalPosition,
2138 if (m_env.numSubEnvironments() > 1) {
2139 V subMeanInter0Mean (m_vectorSpace.zeroVector());
2140 V subMeanInter0Clt95 (m_vectorSpace.zeroVector());
2141 V subMeanInter0Empirical90(m_vectorSpace.zeroVector());
2142 V subMeanInter0Min (m_vectorSpace.zeroVector());
2143 V subMeanInter0Max (m_vectorSpace.zeroVector());
2144 this->subMeanInter0MonitorRun(currentMonitoredFinalPosition,
2147 subMeanInter0Empirical90,
2150 this->subMeanInter0MonitorStore(i,
2151 currentMonitoredFinalPosition,
2154 subMeanInter0Empirical90,
2158 V unifiedMeanVec (m_vectorSpace.zeroVector());
2159 V unifiedMeanCltStd(m_vectorSpace.zeroVector());
2160 this->unifiedMeanMonitorRun(currentMonitoredFinalPosition,
2163 this->unifiedMeanMonitorStore(i,
2164 currentMonitoredFinalPosition,
2171 this->subMeanMonitorWrite(*passedOfs);
2172 if (m_env.numSubEnvironments() > 1) {
2173 this->subMeanInter0MonitorWrite(*passedOfs);
2174 this->unifiedMeanMonitorWrite(*passedOfs);
2178 this->subMeanMonitorFree();
2179 if (m_env.numSubEnvironments() > 1) {
2180 this->subMeanInter0MonitorFree();
2181 this->unifiedMeanMonitorFree();
2185 if (m_env.subDisplayFile()) {
2186 *m_env.subDisplayFile() <<
"Mean evolution took " << tmpRunTime
2193 #endif //#ifdef UQ_CODE_HAS_MONITORS
2200 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2201 template<
class V,
class M>
2203 BaseVectorSequence<V,M>::computeBMM(
2204 const SequenceStatisticalOptions& statisticalOptions,
2205 const std::vector<unsigned int>& initialPosForStatistics,
2206 std::ofstream* passedOfs)
2209 computeMeanStacc(statisticalOptions.m_ov, initialPosForStatistics, passedOfs);
2212 template<
class V,
class M>
2214 BaseVectorSequence<V,M>::computeMeanStacc(
2215 const SsOptionsValues& statisticalOptions,
2216 const std::vector<unsigned int>& initialPosForStatistics,
2217 std::ofstream* passedOfs)
2220 struct timeval timevalTmp;
2221 iRC = gettimeofday(&timevalTmp, NULL);
2222 double tmpRunTime = 0.;
2224 if (m_env.subDisplayFile()) {
2225 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2226 <<
"\nComputing variance of sample mean through BMM"
2230 if (m_env.subDisplayFile()) {
2231 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeBMM(): lengths for batchs in BMM =";
2232 for (
unsigned int i = 0; i < statisticalOptions.bmmLengths().size(); ++i) {
2233 *m_env.subDisplayFile() <<
" " << statisticalOptions.bmmLengths()[i];
2235 *m_env.subDisplayFile() << std::endl;
2238 TwoDArray<V> _2dArrayOfBMM(initialPosForStatistics.size(),statisticalOptions.bmmLengths().size());
2239 for (
unsigned int i = 0; i < _2dArrayOfBMM.numRows(); ++i) {
2240 for (
unsigned int j = 0; j < _2dArrayOfBMM.numCols(); ++j) {
2241 _2dArrayOfBMM.setLocation(i,j,
new V(m_vectorSpace.zeroVector()) );
2244 V bmmVec(m_vectorSpace.zeroVector());
2245 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2246 unsigned int initialPos = initialPosForStatistics[initialPosId];
2247 for (
unsigned int batchLengthId = 0; batchLengthId < statisticalOptions.bmmLengths().size(); batchLengthId++) {
2248 unsigned int batchLength = statisticalOptions.bmmLengths()[batchLengthId];
2249 this->bmm(initialPos,
2252 _2dArrayOfBMM(initialPosId,batchLengthId) = bmmVec;
2256 if (m_env.subDisplayFile()) {
2257 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2258 *m_env.subDisplayFile() <<
"\nEstimated variance of sample mean, through batch means method, for subchain beginning at position " << initialPosForStatistics[initialPosId]
2259 <<
" (each column corresponds to a batch length)"
2265 *m_env.subDisplayFile() << line;
2266 for (
unsigned int batchLengthId = 0; batchLengthId < statisticalOptions.bmmLengths().size(); batchLengthId++) {
2267 sprintf(line,
"%10s%3d",
2269 statisticalOptions.bmmLengths()[batchLengthId]);
2270 *m_env.subDisplayFile() << line;
2273 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2274 sprintf(line,
"\n%9.9s",
2275 m_vectorSpace.localComponentName(i).c_str() );
2276 *m_env.subDisplayFile() << line;
2277 for (
unsigned int batchLengthId = 0; batchLengthId < statisticalOptions.bmmLengths().size(); batchLengthId++) {
2278 sprintf(line,
"%2s%11.4e",
2280 _2dArrayOfBMM(initialPosId,batchLengthId)[i]);
2281 *m_env.subDisplayFile() << line;
2284 *m_env.subDisplayFile() << std::endl;
2289 if (m_env.subDisplayFile()) {
2290 *m_env.subDisplayFile() <<
"Chain BMM took " << tmpRunTime
2298 template<
class V,
class M>
2300 BaseVectorSequence<V,M>::computeFFT(
2301 const SequenceStatisticalOptions& statisticalOptions,
2302 const std::vector<unsigned int>& initialPosForStatistics,
2303 std::ofstream* passedOfs)
2306 computeBMM(statisticalOptions.m_ov, initialPosForStatistics, passedOfs);
2309 template<
class V,
class M>
2311 BaseVectorSequence<V,M>::computeBMM(
2312 const SsOptionsValues& statisticalOptions,
2313 const std::vector<unsigned int>& initialPosForStatistics,
2314 std::ofstream* passedOfs)
2317 struct timeval timevalTmp;
2318 iRC = gettimeofday(&timevalTmp, NULL);
2319 double tmpRunTime = 0.;
2321 if (m_env.subDisplayFile()) {
2322 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2323 <<
"\nComputing FFT of chain on parameter of id = " << statisticalOptions.fftParamId()
2327 std::vector<std::complex<double> > forwardResult(0,std::complex<double>(0.,0.));
2328 std::vector<std::complex<double> > inverseResult(0,std::complex<double>(0.,0.));
2329 Fft<std::complex<double> > fftObj(m_env);
2330 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2331 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2332 this->fftForward(initialPosition,
2333 statisticalOptions.fftSize(),
2334 statisticalOptions.fftParamId(),
2337 if (statisticalOptions.fftWrite() && passedOfs) {
2338 std::ofstream& ofsvar = *passedOfs;
2339 ofsvar << m_name <<
"_fftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << 1
2340 <<
"," << forwardResult.size()
2343 for (
unsigned int j = 0; j < forwardResult.size(); ++j) {
2344 ofsvar << m_name <<
"_fftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << 1
2346 <<
") = " << forwardResult[j].real()
2347 <<
" + i*" << forwardResult[j].imag()
2353 if (statisticalOptions.fftTestInversion()) {
2354 fftObj.inverse(forwardResult,
2355 statisticalOptions.fftSize(),
2357 if (statisticalOptions.fftWrite() && passedOfs) {
2358 std::ofstream& ofsvar = *passedOfs;
2359 ofsvar << m_name <<
"_iftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << 1
2360 <<
"," << inverseResult.size()
2363 for (
unsigned int j = 0; j < inverseResult.size(); ++j) {
2364 ofsvar << m_name <<
"_iftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << 1
2366 <<
") = " << inverseResult[j].real()
2367 <<
" + i*" << inverseResult[j].imag()
2376 if (m_env.subDisplayFile()) {
2377 *m_env.subDisplayFile() <<
"Chain FFT took " << tmpRunTime
2385 template<
class V,
class M>
2387 BaseVectorSequence<V,M>::computePSD(
2388 const SequenceStatisticalOptions& statisticalOptions,
2389 const std::vector<unsigned int>& initialPosForStatistics,
2390 std::ofstream* passedOfs)
2393 computeFFT(statisticalOptions.m_ov, initialPosForStatistics, passedOfs);
2396 template<
class V,
class M>
2398 BaseVectorSequence<V,M>::computeFFT(
2399 const SsOptionsValues& statisticalOptions,
2400 const std::vector<unsigned int>& initialPosForStatistics,
2401 std::ofstream* passedOfs)
2404 struct timeval timevalTmp;
2405 iRC = gettimeofday(&timevalTmp, NULL);
2406 double tmpRunTime = 0.;
2408 if (m_env.subDisplayFile()) {
2409 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2410 <<
"\nComputing PSD of chain on parameter of id = " << statisticalOptions.psdParamId()
2414 std::vector<double> psdResult(0,0.);
2415 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2416 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2417 this->psd(initialPosition,
2418 statisticalOptions.psdNumBlocks(),
2419 statisticalOptions.psdHopSizeRatio(),
2420 statisticalOptions.psdParamId(),
2423 if (statisticalOptions.psdWrite() && passedOfs) {
2424 std::ofstream& ofsvar = *passedOfs;
2425 ofsvar << m_name <<
"_psdInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << 1
2426 <<
"," << psdResult.size()
2429 for (
unsigned int j = 0; j < psdResult.size(); ++j) {
2430 ofsvar << m_name <<
"_psdInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << 1
2432 <<
") = " << psdResult[j]
2440 if (m_env.subDisplayFile()) {
2441 *m_env.subDisplayFile() <<
"Chain PSD took " << tmpRunTime
2449 template<
class V,
class M>
2451 BaseVectorSequence<V,M>::computePSDAtZero(
2452 const SequenceStatisticalOptions& statisticalOptions,
2453 const std::vector<unsigned int>& initialPosForStatistics,
2454 std::ofstream* passedOfs)
2457 computePSD(statisticalOptions.m_ov, initialPosForStatistics, passedOfs);
2460 template<
class V,
class M>
2462 BaseVectorSequence<V,M>::computePSD(
2463 const SsOptionsValues& statisticalOptions,
2464 const std::vector<unsigned int>& initialPosForStatistics,
2465 std::ofstream* passedOfs)
2468 struct timeval timevalTmp;
2469 iRC = gettimeofday(&timevalTmp, NULL);
2470 double tmpRunTime = 0.;
2472 if (m_env.subDisplayFile()) {
2473 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2474 <<
"\nComputing PSD at frequency zero for all parameters"
2478 TwoDArray<V> _2dArrayOfPSDAtZero(initialPosForStatistics.size(),statisticalOptions.psdAtZeroNumBlocks().size());
2479 for (
unsigned int i = 0; i < _2dArrayOfPSDAtZero.numRows(); ++i) {
2480 for (
unsigned int j = 0; j < _2dArrayOfPSDAtZero.numCols(); ++j) {
2481 _2dArrayOfPSDAtZero.setLocation(i,j,
new V(m_vectorSpace.zeroVector()) );
2484 V psdVec(m_vectorSpace.zeroVector());
2485 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2486 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2487 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2488 unsigned int numBlocks = statisticalOptions.psdAtZeroNumBlocks()[numBlocksId];
2489 this->psdAtZero(initialPosition,
2491 statisticalOptions.psdAtZeroHopSizeRatio(),
2493 _2dArrayOfPSDAtZero(initialPosId,numBlocksId) = psdVec;
2498 if ((statisticalOptions.psdAtZeroDisplay()) && (m_env.subDisplayFile())) {
2499 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2500 unsigned int initialPos = initialPosForStatistics[initialPosId];
2501 *m_env.subDisplayFile() <<
"\nComputed PSD at frequency zero for subchain beginning at position " << initialPos
2502 <<
", so effective data size = " << this->subSequenceSize() - initialPos
2503 <<
" (each column corresponds to a number of blocks)"
2509 *m_env.subDisplayFile() << line;
2510 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2511 sprintf(line,
"%10s%3d",
2513 statisticalOptions.psdAtZeroNumBlocks()[numBlocksId]);
2514 *m_env.subDisplayFile() << line;
2517 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2518 sprintf(line,
"\n%9.9s",
2519 m_vectorSpace.localComponentName(i).c_str() );
2520 *m_env.subDisplayFile() << line;
2521 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2522 sprintf(line,
"%2s%11.4e",
2524 _2dArrayOfPSDAtZero(initialPosId,numBlocksId)[i]);
2525 *m_env.subDisplayFile() << line;
2528 *m_env.subDisplayFile() << std::endl;
2533 if ( (m_env.subDisplayFile())) {
2534 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2535 unsigned int initialPos = initialPosForStatistics[initialPosId];
2536 *m_env.subDisplayFile() <<
"\nEstimated variance of sample mean, through psd, for subchain beginning at position " << initialPos
2537 <<
", so effective data size = " << this->subSequenceSize() - initialPos
2538 <<
" (each column corresponds to a number of blocks)"
2544 *m_env.subDisplayFile() << line;
2545 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2546 sprintf(line,
"%10s%3d",
2548 statisticalOptions.psdAtZeroNumBlocks()[numBlocksId]);
2549 *m_env.subDisplayFile() << line;
2552 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2553 sprintf(line,
"\n%9.9s",
2554 m_vectorSpace.localComponentName(i).c_str() );
2555 *m_env.subDisplayFile() << line;
2556 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2557 sprintf(line,
"%2s%11.4e",
2559 2.*M_PI*_2dArrayOfPSDAtZero(initialPosId,numBlocksId)[i]/(
double) (this->subSequenceSize() - initialPos));
2560 *m_env.subDisplayFile() << line;
2563 *m_env.subDisplayFile() << std::endl;
2568 if (m_env.subDisplayFile()) {
2569 *m_env.subDisplayFile() <<
"Chain PSD at frequency zero took " << tmpRunTime
2575 if (statisticalOptions.psdAtZeroWrite() && passedOfs) {
2576 std::ofstream& ofsvar = *passedOfs;
2577 ofsvar << m_name <<
"_psdAtZeroNumBlocks_sub" << m_env.subIdString() <<
" = zeros(" << 1
2578 <<
"," << statisticalOptions.psdAtZeroNumBlocks().size()
2581 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2582 ofsvar << m_name <<
"_psdAtZeroNumBlocks_sub" << m_env.subIdString() <<
"(" << 1
2583 <<
"," << numBlocksId+1
2584 <<
") = " << statisticalOptions.psdAtZeroNumBlocks()[numBlocksId]
2589 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2590 ofsvar << m_name <<
"_psdAtZeroInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << this->vectorSizeLocal()
2591 <<
"," << statisticalOptions.psdAtZeroNumBlocks().size()
2594 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2595 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2596 ofsvar << m_name <<
"_psdAtZeroInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << i+1
2597 <<
"," << numBlocksId+1
2598 <<
") = " << _2dArrayOfPSDAtZero(initialPosId,numBlocksId)[i]
2609 template<
class V,
class M>
2611 BaseVectorSequence<V,M>::computeGeweke(
2612 const SequenceStatisticalOptions& statisticalOptions,
2613 const std::vector<unsigned int>& initialPosForStatistics,
2614 std::ofstream* passedOfs)
2617 computePSDAtZero(statisticalOptions.m_ov, initialPosForStatistics, passedOfs);
2620 template<
class V,
class M>
2622 BaseVectorSequence<V,M>::computePSDAtZero(
2623 const SsOptionsValues& statisticalOptions,
2624 const std::vector<unsigned int>& initialPosForStatistics,
2625 std::ofstream* passedOfs)
2628 struct timeval timevalTmp;
2629 iRC = gettimeofday(&timevalTmp, NULL);
2630 double tmpRunTime = 0.;
2632 if (m_env.subDisplayFile()) {
2633 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2634 <<
"\nComputing Geweke coefficients"
2638 std::vector<V*> vectorOfGeweke(initialPosForStatistics.size(),NULL);
2639 V gewVec(m_vectorSpace.zeroVector());
2640 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2641 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2642 this->geweke(initialPosition,
2643 statisticalOptions.gewekeNaRatio(),
2644 statisticalOptions.gewekeNbRatio(),
2646 vectorOfGeweke[initialPosId] =
new V(gewVec);
2649 if (m_env.subDisplayFile()) {
2650 *m_env.subDisplayFile() <<
"\nComputed Geweke coefficients with 10% and 50% percentages"
2651 <<
" (each column corresponds to a different initial position on the full chain)"
2657 *m_env.subDisplayFile() << line;
2658 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2659 sprintf(line,
"%10s%3d",
2661 initialPosForStatistics[initialPosId]);
2662 *m_env.subDisplayFile() << line;
2665 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2666 sprintf(line,
"\n%9.9s",
2667 m_vectorSpace.localComponentName(i).c_str() );
2668 *m_env.subDisplayFile() << line;
2669 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2670 sprintf(line,
"%2s%11.4e",
2672 (*(vectorOfGeweke[initialPosId]))[i]);
2673 *m_env.subDisplayFile() << line;
2676 *m_env.subDisplayFile() << std::endl;
2680 if (m_env.subDisplayFile()) {
2681 *m_env.subDisplayFile() <<
"Chain Geweke took " << tmpRunTime
2689 template<
class V,
class M>
2691 BaseVectorSequence<V,M>::computeMeanStacc(
2692 const SequenceStatisticalOptions& statisticalOptions,
2693 const std::vector<unsigned int>& initialPosForStatistics,
2694 std::ofstream* passedOfs)
2697 computeGeweke(statisticalOptions.m_ov, initialPosForStatistics, passedOfs);
2700 template<
class V,
class M>
2702 BaseVectorSequence<V,M>::computeGeweke(
2703 const SsOptionsValues& statisticalOptions,
2704 const std::vector<unsigned int>& initialPosForStatistics,
2705 std::ofstream* passedOfs)
2708 struct timeval timevalTmp;
2709 iRC = gettimeofday(&timevalTmp, NULL);
2710 double tmpRunTime = 0.;
2712 if (m_env.subDisplayFile()) {
2713 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2714 <<
"\nComputing mean statistical accuracy"
2718 std::vector<V*> vectorOfMeanStacc(initialPosForStatistics.size(),NULL);
2719 V meanStaccVec(m_vectorSpace.zeroVector());
2720 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2721 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2722 this->meanStacc(initialPosition,
2724 vectorOfMeanStacc[initialPosId] =
new V(meanStaccVec);
2727 if (m_env.subDisplayFile()) {
2728 *m_env.subDisplayFile() <<
"\nComputed mean statistical accuracy"
2729 <<
" (each column corresponds to a different initial position on the full chain)"
2735 *m_env.subDisplayFile() << line;
2736 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2737 sprintf(line,
"%10s%3d",
2739 initialPosForStatistics[initialPosId]);
2740 *m_env.subDisplayFile() << line;
2743 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2744 sprintf(line,
"\n%9.9s",
2745 m_vectorSpace.localComponentName(i).c_str() );
2746 *m_env.subDisplayFile() << line;
2747 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2748 sprintf(line,
"%2s%11.4e",
2750 (*(vectorOfMeanStacc[initialPosId]))[i]);
2751 *m_env.subDisplayFile() << line;
2754 *m_env.subDisplayFile() << std::endl;
2758 if (m_env.subDisplayFile()) {
2759 *m_env.subDisplayFile() <<
"Chain mean statistical accuracy took " << tmpRunTime
2766 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2772 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
2777 unsigned int subNumSamples,
2782 "must provide at least 2 samples to compute correlation matrices");
2795 queso_require_msg(!((numRowsLocal != pqCovMatrix.numRowsLocal()) || (numCols != pqCovMatrix.numCols())),
"inconsistent dimensions for covariance matrix");
2797 queso_require_msg(!((numRowsLocal != pqCorrMatrix.numRowsLocal()) || (numCols != pqCorrMatrix.numCols())),
"inconsistent dimensions for correlation matrix");
2813 for (
unsigned i = 0; i < numRowsLocal; ++i) {
2814 for (
unsigned j = 0; j < numCols; ++j) {
2815 pqCovMatrix(i,j) = 0.;
2818 for (
unsigned k = 0;
k < subNumSamples; ++
k) {
2821 tmpP -= unifiedMeanP;
2824 tmpQ -= unifiedMeanQ;
2826 for (
unsigned i = 0; i < numRowsLocal; ++i) {
2827 for (
unsigned j = 0; j < numCols; ++j) {
2828 pqCovMatrix(i,j) += tmpP[i]*tmpQ[j];
2838 unifiedSampleVarianceP);
2844 unifiedSampleVarianceQ);
2847 double minSampleVarianceP;
2848 double minSampleVarianceQ;
2849 minSampleVarianceP = unifiedSampleVarianceP.getMinValue();
2850 minSampleVarianceQ = unifiedSampleVarianceQ.getMinValue();
2855 if (useOnlyInter0Comm) {
2857 unsigned int unifiedNumSamples = 0;
2859 "ComputeCovCorrMatricesBetweenVectorSequences()",
2860 "failed MPI.Allreduce() for subNumSamples");
2862 for (
unsigned i = 0; i < numRowsLocal; ++i) {
2863 for (
unsigned j = 0; j < numCols; ++j) {
2866 "ComputeCovCorrMatricesBetweenVectorSequences()",
2867 "failed MPI.Allreduce() for a matrix position");
2868 pqCovMatrix(i,j) = aux/((double) (unifiedNumSamples-1));
2872 for (
unsigned i = 0; i < numRowsLocal; ++i) {
2873 for (
unsigned j = 0; j < numCols; ++j) {
2874 pqCorrMatrix(i,j) = pqCovMatrix(i,j)/std::sqrt(unifiedSampleVarianceP[i])/std::sqrt(unifiedSampleVarianceQ[j]);
2875 if (((pqCorrMatrix(i,j) + 1.) < -1.e-8) ||
2876 ((pqCorrMatrix(i,j) - 1.) > 1.e-8)) {
2878 std::cerr <<
"In ComputeCovCorrMatricesBetweenVectorSequences()"
2882 <<
", pqCorrMatrix(i,j)+1 = " << pqCorrMatrix(i,j)+1.
2883 <<
", pqCorrMatrix(i,j)-1 = " << pqCorrMatrix(i,j)-1.
2889 (pqCorrMatrix(i,j), -1. - 1.e-8,
2890 "computed correlation is out of range");
2892 (pqCorrMatrix(i,j), 1. + 1.e-8,
2893 "computed correlation is out of range");
const V & unifiedMaxPlain() const
Finds the maximum value of the unified sequence.
unsigned int vectorSizeLocal() const
Local dimension (size) of the vector space.
const V & unifiedMeanPlain() const
Finds the mean value of the unified sequence.
void append(const BaseVectorSequence< V, M > &src, unsigned int initialPos, unsigned int numPos)
Appends the vector src to this vector.
#define queso_error_msg(msg)
const V & unifiedSampleVariancePlain() const
Finds the variance of a sample of the unified sequence.
double unifiedPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &unifiedPositionsOfMaximum)
Finds the positions where the maximum element occurs in the unified 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...
const BoxSubset< V, M > & unifiedBoxPlain() const
Finds a box subset of the unified-sequence (given by the min and max values of the unified sequence c...
void computeFilterParams(std::ofstream *passedOfs, unsigned int &initialPos, unsigned int &spacing)
Computes the filtering parameters spacing for the sequence of vectors.
unsigned int vectorSizeGlobal() const
Global dimension (size) of the vector space.
void ComputeCovCorrMatricesBetweenVectorSequences(const BaseVectorSequence< P_V, P_M > &subPSeq, const BaseVectorSequence< Q_V, Q_M > &subQSeq, unsigned int subNumSamples, P_M &pqCovMatrix, P_M &pqCorrMatrix)
void deleteStoredVectors()
Deletes all the stored vectors.
double subPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &subPositionsOfMaximum)
Finds the positions where the maximum element occurs in the sub-sequence.
unsigned int dimGlobal() const
Class for a Fast Fourier Transform (FFT) algorithm.
const T & subMaxPlain() const
Finds the maximum value of the sub-sequence of scalars.
const VectorSpace< V, M > & m_vectorSpace
const V & unifiedMedianPlain() const
Finds the median value of the unified sequence.
#define queso_deprecated()
double MiscGetEllapsedSeconds(struct timeval *timeval0)
#define queso_require_less_equal_msg(expr1, expr2, msg)
#define queso_require_msg(asserted, msg)
int inter0Rank() const
Returns the process inter0 rank.
virtual void unifiedMeanExtra(unsigned int initialPos, unsigned int numPos, V &unifiedMeanVec) const =0
Finds the mean value of the unified sequence of numPos positions starting at position initialPos...
const V & subMaxPlain() const
Finds the maximum value of the sub-sequence.
#define queso_require_equal_to_msg(expr1, expr2, msg)
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
virtual void unifiedSampleVarianceExtra(unsigned int initialPos, unsigned int numPos, const V &unifiedMeanVec, V &unifiedSamVec) const =0
Finds the sample variance of the unified sequence, considering numPos positions starting at position ...
#define queso_require_greater_equal_msg(expr1, expr2, msg)
void clear()
Reset the values and the size of the sequence of vectors.
const V & subSampleVariancePlain() const
Finds the variance of a sample of the sub-sequence.
const VectorSpace< V, M > & vectorSpace() const
Vector space; access to protected attribute VectorSpace<V,M>& m_vectorSpace.
Class for matrix operations using GSL library.
const V & unifiedMinPlain() const
Finds the minimum value of the unified sequence.
const V & subMinPlain() const
Finds the minimum value of the sub-sequence.
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
const V & zeroVector() const
Returns a vector filled with zeros.
BaseVectorSequence(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
const std::string & name() const
Access to protected attribute m_name: name of the sequence of vectors.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
virtual void getPositionValues(unsigned int posId, V &vec) const =0
Gets the values of the sequence at position posId and stores them at vec. See template specialization...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
Class representing a subset of a vector space shaped like a hypercube.
unsigned int numOfProcsForStorage() const
Returns total number of processes.
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
const BoxSubset< V, M > & subBoxPlain() const
Finds a box subset of the sub-sequence (given by its min and max values calculated via subMinPlain an...
A class representing a vector space.
virtual ~BaseVectorSequence()
Destructor.
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...
int worldRank() const
Returns the process world rank.
void copy(const BaseVectorSequence< V, M > &src)
Copies vector sequence src to this.
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
unsigned int dimLocal() const
#define queso_require_greater_msg(expr1, expr2, msg)
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
const V & subMeanPlain() const
Finds the mean value of the sub-sequence.
Base class for handling vector and array samples (sequence of vectors or arrays). ...
void MiscComputePositionsBetweenMinMax(V minValues, V maxValues, std::vector< V * > &positions)
unsigned int unifiedSequenceSize() const
Calculates the size of the unified sequence of vectors.
const V & subMedianPlain() const
Finds the median value of the sub-sequence.
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization.