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();
80 "BaseVectorSequence<V,M>::unifiedSequenceSize()",
81 "failed MPI.Allreduce() for unifiedSequenceSize()");
85 unifiedNumSamples = this->subSequenceSize();
91 "BaseVectorSequence<V,M>::unifiedSequenceSize()",
92 "parallel vectors not supported yet");
95 return unifiedNumSamples;
98 template <
class V,
class M>
102 return m_vectorSpace.dimLocal();
105 template <
class V,
class M>
109 return m_vectorSpace.dimGlobal();
112 template <
class V,
class M>
116 return m_vectorSpace;
119 template <
class V,
class M>
126 template <
class V,
class M>
134 template <
class V,
class M>
138 unsigned int numPos = this->subSequenceSize();
140 this->resetValues(0,numPos);
141 this->resizeSequence(0);
147 template <
class V,
class M>
151 if (m_subMinPlain == NULL) {
152 m_subMinPlain = m_vectorSpace.newVector();
153 if (m_subMaxPlain == NULL) m_subMaxPlain = m_vectorSpace.newVector();
154 subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
157 return *m_subMinPlain;
160 template <
class V,
class M>
164 if (m_unifiedMinPlain == NULL) {
165 m_unifiedMinPlain = m_vectorSpace.newVector();
166 if (m_unifiedMaxPlain == NULL) m_unifiedMaxPlain = m_vectorSpace.newVector();
167 unifiedMinMaxExtra(0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
170 return *m_unifiedMinPlain;
173 template <
class V,
class M>
177 if (m_subMaxPlain == NULL) {
178 if (m_subMinPlain == NULL) m_subMinPlain = m_vectorSpace.newVector();
179 m_subMaxPlain = m_vectorSpace.newVector();
180 subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
183 return *m_subMaxPlain;
186 template <
class V,
class M>
190 if (m_unifiedMaxPlain == NULL) {
191 if (m_unifiedMinPlain == NULL) m_unifiedMinPlain = m_vectorSpace.newVector();
192 m_unifiedMaxPlain = m_vectorSpace.newVector();
193 unifiedMinMaxExtra(0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
196 return *m_unifiedMaxPlain;
199 template <
class V,
class M>
203 if (m_subMeanPlain == NULL) {
204 m_subMeanPlain = m_vectorSpace.newVector();
205 subMeanExtra(0,subSequenceSize(),*m_subMeanPlain);
208 return *m_subMeanPlain;
211 template <
class V,
class M>
215 if (m_unifiedMeanPlain == NULL) {
216 m_unifiedMeanPlain = m_vectorSpace.newVector();
217 unifiedMeanExtra(0,subSequenceSize(),*m_unifiedMeanPlain);
220 return *m_unifiedMeanPlain;
223 template <
class V,
class M>
227 if (m_subMedianPlain == NULL) {
228 m_subMedianPlain = m_vectorSpace.newVector();
229 subMedianExtra(0, subSequenceSize(), *m_subMedianPlain);
232 return *m_subMedianPlain;
235 template <
class V,
class M>
239 if (m_unifiedMedianPlain == NULL) {
240 m_unifiedMedianPlain = m_vectorSpace.newVector();
241 unifiedMedianExtra(0, subSequenceSize(), *m_unifiedMedianPlain);
244 return *m_unifiedMedianPlain;
247 template <
class V,
class M>
251 if (m_subSampleVariancePlain == NULL) {
252 m_subSampleVariancePlain = m_vectorSpace.newVector();
253 subSampleVarianceExtra(0,subSequenceSize(),subMeanPlain(),*m_subSampleVariancePlain);
256 return *m_subSampleVariancePlain;
259 template <
class V,
class M>
263 if (m_unifiedSampleVariancePlain == NULL) {
264 m_unifiedSampleVariancePlain = m_vectorSpace.newVector();
265 unifiedSampleVarianceExtra(0,subSequenceSize(),unifiedMeanPlain(),*m_unifiedSampleVariancePlain);
268 return *m_unifiedSampleVariancePlain;
271 template <
class V,
class M>
275 if (m_subBoxPlain == NULL) {
279 this->subMaxPlain());
282 return *m_subBoxPlain;
285 template <
class V,
class M>
289 if (m_unifiedBoxPlain == NULL) {
292 this->unifiedMinPlain(),
293 this->unifiedMaxPlain());
296 return *m_unifiedBoxPlain;
299 template <
class V,
class M>
304 delete m_subMinPlain;
305 m_subMinPlain = NULL;
307 if (m_unifiedMinPlain) {
308 delete m_unifiedMinPlain;
309 m_unifiedMinPlain = NULL;
312 delete m_subMaxPlain;
313 m_subMaxPlain = NULL;
315 if (m_unifiedMaxPlain) {
316 delete m_unifiedMaxPlain;
317 m_unifiedMaxPlain = NULL;
319 if (m_subMeanPlain) {
320 delete m_subMeanPlain;
321 m_subMeanPlain = NULL;
323 if (m_unifiedMeanPlain) {
324 delete m_unifiedMeanPlain;
325 m_unifiedMeanPlain = NULL;
327 if (m_subMedianPlain) {
328 delete m_subMedianPlain;
329 m_subMedianPlain = NULL;
331 if (m_unifiedMedianPlain) {
332 delete m_unifiedMedianPlain;
333 m_unifiedMedianPlain = NULL;
335 if (m_subSampleVariancePlain) {
336 delete m_subSampleVariancePlain;
337 m_subSampleVariancePlain = NULL;
339 if (m_unifiedSampleVariancePlain) {
340 delete m_unifiedSampleVariancePlain;
341 m_unifiedSampleVariancePlain = NULL;
344 delete m_subBoxPlain;
345 m_subBoxPlain = NULL;
347 if (m_unifiedBoxPlain) {
348 delete m_unifiedBoxPlain;
349 m_unifiedBoxPlain = NULL;
355 template <
class V,
class M>
359 unsigned int initialPos,
364 "BaseVectorSequence<V,M>::append()",
365 "initialPos is too big");
369 "BaseVectorSequence<V,M>::append()",
370 "numPos is too big");
372 this->deleteStoredVectors();
373 unsigned int currentSize = this->subSequenceSize();
374 this->resizeSequence(currentSize+numPos);
376 for (
unsigned int i = 0; i < numPos; ++i) {
378 this->setPositionValues(currentSize+i,tmpVec);
384 template <
class V,
class M>
390 if (m_env.subDisplayFile()) {
391 *m_env.subDisplayFile() <<
"Entering BaseVectorSequence<V,M>::subPositionsOfMaximum()"
392 <<
": subCorrespondingScalarValues,subSequenceSize() = " << subCorrespondingScalarValues.
subSequenceSize()
393 <<
", this->subSequenceSize = " << this->subSequenceSize()
399 "BaseVectorSequence<V,M>::subPositionsOfMaximum()",
402 double subMaxValue = subCorrespondingScalarValues.
subMaxPlain();
405 unsigned int subNumPos = 0;
406 for (
unsigned int i = 0; i < iMax; ++i) {
407 if (subCorrespondingScalarValues[i] == subMaxValue) {
412 V tmpVec(this->vectorSpace().zeroVector());
415 for (
unsigned int i = 0; i < iMax; ++i) {
416 if (subCorrespondingScalarValues[i] == subMaxValue) {
417 this->getPositionValues (i,tmpVec);
423 if (m_env.subDisplayFile()) {
424 *m_env.subDisplayFile() <<
"Leaving BaseVectorSequence<V,M>::subPositionsOfMaximum()"
431 template <
class V,
class M>
437 if (m_env.subDisplayFile()) {
438 *m_env.subDisplayFile() <<
"Entering BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()"
439 <<
": subCorrespondingScalarValues,subSequenceSize() = " << subCorrespondingScalarValues.
subSequenceSize()
440 <<
", this->subSequenceSize = " << this->subSequenceSize()
446 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
450 double subMaxValue = subCorrespondingScalarValues.
subMaxPlain();
455 double unifiedMaxValue;
456 std::vector<double> sendbufPos(1,0.);
457 for (
unsigned int i = 0; i < sendbufPos.size(); ++i) {
458 sendbufPos[i] = subMaxValue;
461 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
462 "failed MPI.Allreduce() for max");
468 for (
unsigned int i = 0; i < iMax; ++i) {
469 if (subCorrespondingScalarValues[i] == unifiedMaxValue) {
476 V tmpVec(this->vectorSpace().zeroVector());
479 for (
unsigned int i = 0; i < iMax; ++i) {
481 if (subCorrespondingScalarValues[i] == unifiedMaxValue) {
482 this->getPositionValues (i,tmpVec);
489 std::vector<int> auxBuf(1,0);
490 int unifiedNumPos = 0;
491 auxBuf[0] = subNumPos;
493 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
494 "failed MPI.Allreduce() for sum");
503 unsigned int Np = (
unsigned int) m_env.inter0Comm().NumProc();
505 std::vector<int> recvcntsPos(Np,0);
507 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
508 "failed MPI.Gatherv()");
509 if (m_env.inter0Rank() == 0) {
512 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
513 "failed MPI.Gather() result at proc 0 (recvcntsPos[0])");
517 std::vector<int> displsPos(Np,0);
518 for (
unsigned int nodeId = 1; nodeId < Np; ++nodeId) {
519 displsPos[nodeId] = displsPos[nodeId-1] + recvcntsPos[nodeId-1];
521 if (m_env.inter0Rank() == 0) {
524 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
525 "failed MPI.Gather() result at proc 0 (unifiedNumPos)");
536 unsigned int dimSize = m_vectorSpace.dimLocal();
537 int subNumDbs = subNumPos * dimSize;
538 std::vector<int> recvcntsDbs(Np,0);
540 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
541 "failed MPI.Gatherv()");
542 if (m_env.inter0Rank() == 0) {
545 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
546 "failed MPI.Gather() result at proc 0 (recvcntsDbs[0])");
551 std::vector<int> displsDbs(Np,0);
552 for (
unsigned int nodeId = 1; nodeId < Np; ++nodeId) {
553 displsDbs[nodeId] = displsDbs[nodeId-1] + recvcntsDbs[nodeId-1];
555 if (m_env.inter0Rank() == 0) {
556 UQ_FATAL_TEST_MACRO(((
int) (unifiedNumPos*dimSize)) != (displsDbs[Np - 1] + recvcntsDbs[Np - 1]),
558 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
559 "failed MPI.Gather() result at proc 0 (unifiedNumPos*dimSize)");
567 std::vector<double> sendbufDbs(subNumDbs,0.);
568 for (
unsigned int i = 0; i < (
unsigned int) subNumPos; ++i) {
570 for (
unsigned int j = 0; j < dimSize; ++j) {
571 sendbufDbs[i*dimSize + j] = tmpVec[j];
575 std::vector<double> recvbufDbs(unifiedNumPos * dimSize);
578 m_env.inter0Comm().Gatherv((
void *) &sendbufDbs[0], (
int) subNumDbs,
RawValue_MPI_DOUBLE, (
void *) &recvbufDbs[0], (
int *) &recvcntsDbs[0], (
int *) &displsDbs[0],
RawValue_MPI_DOUBLE, 0,
579 "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
580 "failed MPI.Gatherv()");
587 if (m_env.inter0Rank() == (int) 0) {
588 for (
unsigned int i = 0; i < (
unsigned int) unifiedNumPos; ++i) {
589 for (
unsigned int j = 0; j < dimSize; ++j) {
590 tmpVec[j] = recvbufDbs[i*dimSize + j];
601 if (m_env.subDisplayFile()) {
602 *m_env.subDisplayFile() <<
"Leaving BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()"
606 return unifiedMaxValue;
609 template <
class V,
class M>
613 V gaussianVector(m_vectorSpace.zeroVector());
614 for (
unsigned int j = 0; j < this->subSequenceSize(); ++j) {
615 gaussianVector.cwSetGaussian(meanVec,stdDevVec);
616 this->setPositionValues(j,gaussianVector);
619 this->deleteStoredVectors();
624 template <
class V,
class M>
628 V uniformVector(m_vectorSpace.zeroVector());
629 for (
unsigned int j = 0; j < this->subSequenceSize(); ++j) {
630 uniformVector.cwSetUniform(aVec,bVec);
631 this->setPositionValues(j,uniformVector);
634 this->deleteStoredVectors();
639 template<
class V,
class M>
642 std::ofstream* passedOfs,
643 unsigned int& initialPos,
644 unsigned int& spacing)
646 if (m_env.subDisplayFile()) {
647 *m_env.subDisplayFile() <<
"\n"
648 <<
"\n-----------------------------------------------------"
649 <<
"\n Computing filter parameters for chain '" << m_name <<
"' ..."
650 <<
"\n-----------------------------------------------------"
655 bool okSituation = ((passedOfs == NULL ) ||
656 ((passedOfs != NULL) && (m_env.subRank() >= 0)));
659 "BaseVectorSequence<V,M>::computeFilterParams()",
660 "unexpected combination of file pointer and subRank");
665 if (m_env.subDisplayFile()) {
666 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
667 <<
"\n Finished computing filter parameters for chain '" << m_name <<
"'"
668 <<
": initialPos = " << initialPos
669 <<
", spacing = " << spacing
670 <<
"\n-----------------------------------------------------"
679 template <
class V,
class M>
687 "SequenceOfVectors<V,M>::copy()",
688 "incompatible vector space dimensions");
691 this->deleteStoredVectors();
701 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
702 template<
class V,
class M>
705 const SequenceStatisticalOptions& statisticalOptions,
706 std::ofstream* passedOfs)
708 if (m_env.subDisplayFile()) {
709 *m_env.subDisplayFile() <<
"\n"
710 <<
"\n-----------------------------------------------------"
711 <<
"\n Computing statistics for chain '" << m_name <<
"' ..."
712 <<
"\n-----------------------------------------------------"
717 bool okSituation = ((passedOfs == NULL ) ||
718 ((passedOfs != NULL) && (m_env.subRank() >= 0)));
721 "BaseVectorSequence<V,M>::computeStatistics()",
722 "unexpected combination of file pointer and subRank");
725 struct timeval timevalTmp;
726 iRC = gettimeofday(&timevalTmp, NULL);
727 double tmpRunTime = 0.;
730 std::vector<unsigned int> initialPosForStatistics(statisticalOptions.initialDiscardedPortions().size(),0);
731 for (
unsigned int i = 0; i < initialPosForStatistics.size(); ++i) {
732 initialPosForStatistics[i] = (
unsigned int) (statisticalOptions.initialDiscardedPortions()[i] * (double) (this->subSequenceSize()-1));
733 if (m_env.subDisplayFile()) {
734 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeStatistics()"
735 <<
": statisticalOptions.initialDiscardedPortions()[" << i <<
"] = " << statisticalOptions.initialDiscardedPortions()[i]
736 <<
", initialPosForStatistics[" << i <<
"] = " << initialPosForStatistics[i]
740 if (m_env.subDisplayFile()) {
741 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeStatistics(): initial positions for statistics =";
742 for (
unsigned int i = 0; i < initialPosForStatistics.size(); ++i) {
743 *m_env.subDisplayFile() <<
" " << initialPosForStatistics[i];
745 *m_env.subDisplayFile() << std::endl;
751 this->computeMeanVars(statisticalOptions,
758 #ifdef UQ_CODE_HAS_MONITORS
759 if (statisticalOptions.meanMonitorPeriod() != 0) {
760 this->computeMeanEvolution(statisticalOptions,
768 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
769 if ((statisticalOptions.bmmRun() ) &&
770 (initialPosForStatistics.size() > 0) &&
771 (statisticalOptions.bmmLengths().size() > 0)) {
772 this->computeBMM(statisticalOptions,
773 initialPosForStatistics,
780 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
781 if ((statisticalOptions.fftCompute() ) &&
782 (initialPosForStatistics.size() > 0)) {
783 this->computeFFT(statisticalOptions,
784 initialPosForStatistics,
791 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
792 if ((statisticalOptions.psdCompute() ) &&
793 (initialPosForStatistics.size() > 0)) {
794 this->computePSD(statisticalOptions,
795 initialPosForStatistics,
802 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
803 if ((statisticalOptions.psdAtZeroCompute() ) &&
804 (initialPosForStatistics.size() > 0) &&
805 (statisticalOptions.psdAtZeroNumBlocks().size() > 0)) {
806 this->computePSDAtZero(statisticalOptions,
807 initialPosForStatistics,
814 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
815 if ((statisticalOptions.gewekeCompute()) &&
816 (initialPosForStatistics.size() > 0)) {
817 this->computeGeweke(statisticalOptions,
818 initialPosForStatistics,
825 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
826 if ((statisticalOptions.meanStaccCompute()) &&
827 (initialPosForStatistics.size() > 0 )) {
828 this->computeMeanStacc(statisticalOptions,
829 initialPosForStatistics,
834 std::vector<unsigned int> lagsForCorrs(statisticalOptions.autoCorrNumLags(),1);
835 for (
unsigned int i = 1; i < lagsForCorrs.size(); ++i) {
836 lagsForCorrs[i] = statisticalOptions.autoCorrSecondLag() + (i-1)*statisticalOptions.autoCorrLagSpacing();
842 if ((statisticalOptions.autoCorrComputeViaDef()) &&
843 (initialPosForStatistics.size() > 0 ) &&
844 (lagsForCorrs.size() > 0 )) {
845 this->computeAutoCorrViaDef(statisticalOptions,
846 initialPosForStatistics,
854 if ((statisticalOptions.autoCorrComputeViaFft()) &&
855 (initialPosForStatistics.size() > 0 ) &&
856 (lagsForCorrs.size() > 0 )) {
857 this->computeAutoCorrViaFFT(statisticalOptions,
858 initialPosForStatistics,
866 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
867 if ((statisticalOptions.histCompute ()) ||
868 (statisticalOptions.cdfStaccCompute()) ||
869 (statisticalOptions.kdeCompute ())) {
871 if (statisticalOptions.kdeCompute())
873 this->computeHistCdfstaccKde(statisticalOptions,
880 if ((statisticalOptions.covMatrixCompute ()) ||
881 (statisticalOptions.corrMatrixCompute())) {
882 this->computeCovCorrMatrices(statisticalOptions,
887 if (m_env.subDisplayFile()) {
888 *m_env.subDisplayFile() <<
"All statistics of chain '" << m_name <<
"'"
889 <<
" took " << tmpRunTime
894 if (m_env.subDisplayFile()) {
895 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
896 <<
"\n Finished computing statistics for chain '" << m_name <<
"'"
897 <<
"\n-----------------------------------------------------"
905 template<
class V,
class M>
907 BaseVectorSequence<V,M>::computeMeanVars(
908 const SequenceStatisticalOptions& statisticalOptions,
909 std::ofstream* passedOfs,
916 struct timeval timevalTmp;
917 iRC = gettimeofday(&timevalTmp, NULL);
918 double tmpRunTime = 0.;
920 if (m_env.subDisplayFile()) {
921 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
922 <<
"\nComputing mean, sample variance and population variance"
926 V subChainMean(m_vectorSpace.zeroVector());
927 this->subMeanExtra(0,
928 this->subSequenceSize(),
931 V subChainMedian(m_vectorSpace.zeroVector());
932 this->subMedianExtra(0,
933 this->subSequenceSize(),
936 V subChainSampleVariance(m_vectorSpace.zeroVector());
937 this->subSampleVarianceExtra(0,
938 this->subSequenceSize(),
940 subChainSampleVariance);
942 if ((m_env.displayVerbosity() >= 5) && (m_env.subDisplayFile())) {
943 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeMeanVars()"
944 <<
": subChainMean.sizeLocal() = " << subChainMean.sizeLocal()
945 <<
", subChainMean = " << subChainMean
946 <<
", subChainMedian = " << subChainMedian
947 <<
", subChainSampleVariance.sizeLocal() = " << subChainSampleVariance.sizeLocal()
948 <<
", subChainSampleVariance = " << subChainSampleVariance
952 V estimatedVarianceOfSampleMean(subChainSampleVariance);
953 estimatedVarianceOfSampleMean /= (double) this->subSequenceSize();
954 bool savedVectorPrintState = estimatedVarianceOfSampleMean.getPrintHorizontally();
955 estimatedVarianceOfSampleMean.setPrintHorizontally(
false);
956 if (m_env.subDisplayFile()) {
957 *m_env.subDisplayFile() <<
"\nEstimated variance of sample mean for the whole chain '" << m_name <<
"'"
958 <<
", under independence assumption:"
960 << estimatedVarianceOfSampleMean
963 estimatedVarianceOfSampleMean.setPrintHorizontally(savedVectorPrintState);
965 V estimatedStdOfSampleMean(estimatedVarianceOfSampleMean);
966 estimatedStdOfSampleMean.cwSqrt();
967 savedVectorPrintState = estimatedStdOfSampleMean.getPrintHorizontally();
968 estimatedStdOfSampleMean.setPrintHorizontally(
false);
969 if (m_env.subDisplayFile()) {
970 *m_env.subDisplayFile() <<
"\nEstimated standard deviation of sample mean for the whole chain '" << m_name <<
"'"
971 <<
", under independence assumption:"
973 << estimatedStdOfSampleMean
976 estimatedStdOfSampleMean.setPrintHorizontally(savedVectorPrintState);
978 V subChainPopulationVariance(m_vectorSpace.zeroVector());
979 this->subPopulationVariance(0,
980 this->subSequenceSize(),
982 subChainPopulationVariance);
985 if (m_env.subDisplayFile()) {
986 *m_env.subDisplayFile() <<
"Sub Mean, median, and variances took " << tmpRunTime
991 if (m_env.subDisplayFile()) {
992 *m_env.subDisplayFile() <<
"\nSub mean, median, sample std, population std"
995 sprintf(line,
"%s%4s%s%6s%s%9s%s%9s%s",
1005 *m_env.subDisplayFile() << line;
1007 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1008 sprintf(line,
"\n%8.8s%2s%11.4e%2s%11.4e%2s%11.4e%2s%11.4e",
1009 m_vectorSpace.localComponentName(i).c_str(),
1015 std::sqrt(subChainSampleVariance[i]),
1017 std::sqrt(subChainPopulationVariance[i]));
1018 *m_env.subDisplayFile() << line;
1020 *m_env.subDisplayFile() << std::endl;
1023 if (subMeanPtr ) *subMeanPtr = subChainMean;
1024 if (subMedianPtr ) *subMedianPtr = subChainMedian;
1025 if (subSampleVarPtr) *subSampleVarPtr = subChainSampleVariance;
1026 if (subPopulVarPtr ) *subPopulVarPtr = subChainPopulationVariance;
1028 if (m_env.numSubEnvironments() > 1) {
1030 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1031 V unifiedChainMean(m_vectorSpace.zeroVector());
1032 this->unifiedMeanExtra(0,
1033 this->subSequenceSize(),
1036 V unifiedChainMedian(m_vectorSpace.zeroVector());
1037 this->unifiedMedianExtra(0,
1038 this->subSequenceSize(),
1039 unifiedChainMedian);
1041 V unifiedChainSampleVariance(m_vectorSpace.zeroVector());
1042 this->unifiedSampleVarianceExtra(0,
1043 this->subSequenceSize(),
1045 unifiedChainSampleVariance);
1047 V unifiedChainPopulationVariance(m_vectorSpace.zeroVector());
1048 this->unifiedPopulationVariance(0,
1049 this->subSequenceSize(),
1051 unifiedChainPopulationVariance);
1053 if (m_env.inter0Rank() == 0) {
1054 if (m_env.subDisplayFile()) {
1055 *m_env.subDisplayFile() <<
"\nUnif mean, median, sample std, population std"
1058 sprintf(line,
"%s%4s%s%6s%s%9s%s%9s%s",
1068 *m_env.subDisplayFile() << line;
1070 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1071 sprintf(line,
"\n%8.8s%2s%11.4e%2s%11.4e%2s%11.4e%2s%11.4e",
1072 m_vectorSpace.localComponentName(i).c_str(),
1074 unifiedChainMean[i],
1076 unifiedChainMedian[i],
1078 std::sqrt(unifiedChainSampleVariance[i]),
1080 std::sqrt(unifiedChainPopulationVariance[i]));
1081 *m_env.subDisplayFile() << line;
1083 *m_env.subDisplayFile() << std::endl;
1090 "BaseVectorSequence<V,M>::computeMeanVars()",
1091 "unified min-max writing, parallel vectors not supported yet");
1095 if (m_env.subDisplayFile()) {
1096 *m_env.subDisplayFile() <<
"\nEnded computing mean, sample variance and population variance"
1097 <<
"\n-----------------------------------------------------"
1104 template<
class V,
class M>
1106 BaseVectorSequence<V,M>::computeAutoCorrViaDef(
1107 const SequenceStatisticalOptions& 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 struct timeval timevalTmp;
1236 iRC = gettimeofday(&timevalTmp, NULL);
1237 double tmpRunTime = 0.;
1239 if (m_env.subDisplayFile()) {
1240 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1241 <<
"\nComputing autocorrelation coefficients (via fft)"
1245 if (statisticalOptions.autoCorrDisplay() && (m_env.subDisplayFile())) {
1246 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeAutoCorrViaFFT(): lags for autocorrelation (via fft) = ";
1247 for (
unsigned int i = 0; i < lagsForCorrs.size(); ++i) {
1248 *m_env.subDisplayFile() <<
" " << lagsForCorrs[i];
1250 *m_env.subDisplayFile() << std::endl;
1253 TwoDArray<V> _2dArrayOfAutoCorrs(initialPosForStatistics.size(),lagsForCorrs.size());
1254 for (
unsigned int i = 0; i < _2dArrayOfAutoCorrs.numRows(); ++i) {
1255 for (
unsigned int j = 0; j < _2dArrayOfAutoCorrs.numCols(); ++j) {
1256 _2dArrayOfAutoCorrs.setLocation(i,j,
new V(m_vectorSpace.zeroVector()) );
1259 std::vector<V*> corrVecs(lagsForCorrs.size(),NULL);
1260 std::vector<V*> corrSumVecs(initialPosForStatistics.size(),NULL);
1261 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1262 corrSumVecs[initialPosId] =
new V(m_vectorSpace.zeroVector()) ;
1263 unsigned int initialPos = initialPosForStatistics[initialPosId];
1264 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1265 corrVecs[lagId] =
new V(m_vectorSpace.zeroVector()) ;
1267 if (m_env.subDisplayFile()) {
1268 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeAutoCorrViaFFT()"
1269 <<
": about to call chain.autoCorrViaFft()"
1270 <<
" with initialPos = " << initialPos
1271 <<
", numPos = " << this->subSequenceSize()-initialPos
1272 <<
", lagsForCorrs.size() = " << lagsForCorrs.size()
1273 <<
", corrVecs.size() = " << corrVecs.size()
1276 this->autoCorrViaFft(initialPos,
1277 this->subSequenceSize()-initialPos,
1280 this->autoCorrViaFft(initialPos,
1281 this->subSequenceSize()-initialPos,
1282 (
unsigned int) (1.0 * (
double) (this->subSequenceSize()-initialPos)),
1283 *corrSumVecs[initialPosId]);
1284 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1285 _2dArrayOfAutoCorrs(initialPosId,lagId) = *(corrVecs[lagId]);
1288 for (
unsigned int j = 0; j < corrVecs.size(); ++j) {
1289 if (corrVecs[j] != NULL)
delete corrVecs[j];
1292 if (statisticalOptions.autoCorrDisplay()) {
1293 V subChainMean (m_vectorSpace.zeroVector());
1294 V subChainSampleVariance (m_vectorSpace.zeroVector());
1295 V estimatedVarianceOfSampleMean(m_vectorSpace.zeroVector());
1296 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1297 unsigned int initialPos = initialPosForStatistics[initialPosId];
1299 this->subMeanExtra(initialPos,
1300 this->subSequenceSize()-initialPos,
1303 this->subSampleVarianceExtra(initialPos,
1304 this->subSequenceSize()-initialPos,
1306 subChainSampleVariance);
1308 if (m_env.subDisplayFile()) {
1309 *m_env.subDisplayFile() <<
"\nEstimated variance of sample mean, through autocorrelation (via fft), for subchain beginning at position " << initialPosForStatistics[initialPosId]
1312 estimatedVarianceOfSampleMean.cwSet(-1.);
1313 estimatedVarianceOfSampleMean += 2.* (*corrSumVecs[initialPosId]);
1314 estimatedVarianceOfSampleMean *= subChainSampleVariance;
1315 estimatedVarianceOfSampleMean /= (double) (this->subSequenceSize() - initialPos);
1316 bool savedVectorPrintState = estimatedVarianceOfSampleMean.getPrintHorizontally();
1317 estimatedVarianceOfSampleMean.setPrintHorizontally(
false);
1318 if (m_env.subDisplayFile()) {
1319 *m_env.subDisplayFile() << estimatedVarianceOfSampleMean
1322 estimatedVarianceOfSampleMean.setPrintHorizontally(savedVectorPrintState);
1324 if (m_env.subDisplayFile()) {
1325 *m_env.subDisplayFile() <<
"\nComputed autocorrelation coefficients (via fft), for subchain beginning at position " << initialPosForStatistics[initialPosId]
1326 <<
" (each column corresponds to a different lag)"
1332 *m_env.subDisplayFile() << line;
1333 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1334 sprintf(line,
"%10s%3d",
1336 lagsForCorrs[lagId]);
1337 *m_env.subDisplayFile() << line;
1340 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1341 sprintf(line,
"\n%9.9s",
1342 m_vectorSpace.localComponentName(i).c_str() );
1343 *m_env.subDisplayFile() << line;
1344 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1345 sprintf(line,
"%2s%11.4e",
1347 _2dArrayOfAutoCorrs(initialPosId,lagId)[i]);
1348 *m_env.subDisplayFile() << line;
1351 *m_env.subDisplayFile() << std::endl;
1357 if (m_env.subDisplayFile()) {
1358 *m_env.subDisplayFile() <<
"Chain autocorrelation (via fft) took " << tmpRunTime
1364 if (statisticalOptions.autoCorrWrite() && passedOfs) {
1365 std::ofstream& ofsvar = *passedOfs;
1366 ofsvar << m_name <<
"_corrViaFftLags_sub" << m_env.subIdString() <<
" = zeros(" << 1
1367 <<
"," << lagsForCorrs.size()
1370 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1371 ofsvar << m_name <<
"_corrViaFftLags_sub" << m_env.subIdString() <<
"(" << 1
1373 <<
") = " << lagsForCorrs[lagId]
1378 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1379 ofsvar << m_name <<
"_corrViaFftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << this->vectorSizeLocal()
1380 <<
"," << lagsForCorrs.size()
1383 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1384 for (
unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1385 ofsvar << m_name <<
"_corrViaFftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << i+1
1387 <<
") = " << _2dArrayOfAutoCorrs(initialPosId,lagId)[i]
1398 template<
class V,
class M>
1400 BaseVectorSequence<V,M>::computeHistCdfstaccKde(
1401 const SequenceStatisticalOptions& statisticalOptions,
1402 std::ofstream* passedOfs)
1404 if (m_env.subDisplayFile()) {
1405 *m_env.subDisplayFile() <<
"\n"
1406 <<
"\n-----------------------------------------------------"
1407 <<
"\n Computing histogram and/or cdf stacc and/or Kde for chain '" << m_name <<
"' ..."
1408 <<
"\n-----------------------------------------------------"
1414 struct timeval timevalTmp;
1419 double tmpRunTime = 0.;
1420 iRC = gettimeofday(&timevalTmp, NULL);
1421 if (m_env.subDisplayFile()) {
1422 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1423 <<
"\nComputing min and max for histograms and Kde"
1427 V statsMinPositions(m_vectorSpace.zeroVector());
1428 V statsMaxPositions(m_vectorSpace.zeroVector());
1429 this->subMinMaxExtra(0,
1430 this->subSequenceSize(),
1434 if (m_env.subDisplayFile()) {
1435 *m_env.subDisplayFile() <<
"\nComputed min values and max values for chain '" << m_name <<
"'"
1441 *m_env.subDisplayFile() << line;
1443 sprintf(line,
"%9s%s%9s%s",
1448 *m_env.subDisplayFile() << line;
1450 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1451 sprintf(line,
"\n%8.8s",
1452 m_vectorSpace.localComponentName(i).c_str() );
1453 *m_env.subDisplayFile() << line;
1455 sprintf(line,
"%2s%11.4e%2s%11.4e",
1457 statsMinPositions[i],
1459 statsMaxPositions[i]);
1460 *m_env.subDisplayFile() << line;
1462 *m_env.subDisplayFile() << std::endl;
1465 V unifiedStatsMinPositions(statsMinPositions);
1466 V unifiedStatsMaxPositions(statsMaxPositions);
1467 if (m_env.numSubEnvironments() > 1) {
1469 this->unifiedMinMaxExtra(0,
1470 this->subSequenceSize(),
1471 unifiedStatsMinPositions,
1472 unifiedStatsMaxPositions);
1475 if (m_env.subDisplayFile()) {
1476 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1477 if (m_env.inter0Rank() == 0) {
1478 *m_env.subDisplayFile() <<
"\nComputed unified min values and max values for chain '" << m_name <<
"'"
1484 *m_env.subDisplayFile() << line;
1486 sprintf(line,
"%9s%s%9s%s",
1491 *m_env.subDisplayFile() << line;
1493 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1494 sprintf(line,
"\n%8.8s",
1495 m_vectorSpace.localComponentName(i).c_str() );
1496 *m_env.subDisplayFile() << line;
1498 sprintf(line,
"%2s%11.4e%2s%11.4e",
1500 unifiedStatsMinPositions[i],
1502 unifiedStatsMaxPositions[i]);
1503 *m_env.subDisplayFile() << line;
1505 *m_env.subDisplayFile() << std::endl;
1511 "BaseVectorSequence<V,M>::computeHistCdfstaccKde()",
1512 "unified min-max writing, parallel vectors not supported yet");
1519 if (m_env.subDisplayFile()) {
1520 *m_env.subDisplayFile() <<
"Chain min and max took " << tmpRunTime
1528 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
1529 if ((statisticalOptions.histCompute() ) &&
1530 (statisticalOptions.histNumInternalBins() > 0)) {
1532 iRC = gettimeofday(&timevalTmp, NULL);
1533 if (m_env.subDisplayFile()) {
1534 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1535 <<
"\nComputing histograms"
1539 std::string subCoreName_HistCenters((std::string)(
"_HistCenters_sub")+m_env.subIdString());
1540 std::string uniCoreName_HistCenters((std::string)(
"_unifHistCenters_sub")+m_env.subIdString());
1541 if (m_env.numSubEnvironments() == 1) subCoreName_HistCenters = uniCoreName_HistCenters;
1543 std::string subCoreName_HistQuantts((std::string)(
"_HistQuantts_sub")+m_env.subIdString());
1544 std::string uniCoreName_HistQuantts((std::string)(
"_unifHistQuantts_sub")+m_env.subIdString());
1545 if (m_env.numSubEnvironments() == 1) subCoreName_HistQuantts = uniCoreName_HistQuantts;
1547 for (
unsigned int i = 0; i < statsMaxPositions.sizeLocal(); ++i) {
1548 statsMaxPositions[i] *= (1. + 1.e-15);
1552 std::vector<V*> histCentersForAllBins(statisticalOptions.histNumInternalBins()+2,NULL);
1553 std::vector<V*> histQuanttsForAllBins(statisticalOptions.histNumInternalBins()+2,NULL);
1554 this->subHistogram(0,
1557 histCentersForAllBins,
1558 histQuanttsForAllBins);
1562 std::ofstream& ofsvar = *passedOfs;
1563 ofsvar << m_name << subCoreName_HistCenters <<
" = zeros(" << this->vectorSizeLocal()
1564 <<
"," << histCentersForAllBins.size()
1567 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1568 for (
unsigned int j = 0; j < histCentersForAllBins.size(); ++j) {
1569 ofsvar << m_name << subCoreName_HistCenters <<
"(" << i+1
1571 <<
") = " << (*(histCentersForAllBins[j]))[i]
1577 ofsvar << m_name << subCoreName_HistQuantts <<
" = zeros(" << this->vectorSizeLocal()
1578 <<
"," << histQuanttsForAllBins.size()
1581 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1582 for (
unsigned int j = 0; j < histQuanttsForAllBins.size(); ++j) {
1583 ofsvar << m_name << subCoreName_HistQuantts <<
"(" << i+1
1585 <<
") = " << (*(histQuanttsForAllBins[j]))[i]
1592 for (
unsigned int i = 0; i < histQuanttsForAllBins.size(); ++i) {
1593 if (histQuanttsForAllBins[i] != NULL)
delete histQuanttsForAllBins[i];
1595 for (
unsigned int i = 0; i < histCentersForAllBins.size(); ++i) {
1596 if (histCentersForAllBins[i] != NULL)
delete histCentersForAllBins[i];
1599 std::vector<V*> unifiedHistCentersForAllBins(statisticalOptions.histNumInternalBins()+2,NULL);
1600 std::vector<V*> unifiedHistQuanttsForAllBins(statisticalOptions.histNumInternalBins()+2,NULL);
1601 if (m_env.numSubEnvironments() > 1) {
1603 this->unifiedHistogram(0,
1604 unifiedStatsMinPositions,
1605 unifiedStatsMaxPositions,
1606 unifiedHistCentersForAllBins,
1607 unifiedHistQuanttsForAllBins);
1611 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1612 if (m_env.inter0Rank() == 0) {
1613 std::ofstream& ofsvar = *passedOfs;
1614 ofsvar << m_name << uniCoreName_HistCenters <<
" = zeros(" << this->vectorSizeLocal()
1615 <<
"," << unifiedHistCentersForAllBins.size()
1618 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1619 for (
unsigned int j = 0; j < unifiedHistCentersForAllBins.size(); ++j) {
1620 ofsvar << m_name << uniCoreName_HistCenters <<
"(" << i+1
1622 <<
") = " << (*(unifiedHistCentersForAllBins[j]))[i]
1628 ofsvar << m_name << uniCoreName_HistQuantts <<
" = zeros(" << this->vectorSizeLocal()
1629 <<
"," << unifiedHistQuanttsForAllBins.size()
1632 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1633 for (
unsigned int j = 0; j < unifiedHistQuanttsForAllBins.size(); ++j) {
1634 ofsvar << m_name << uniCoreName_HistQuantts <<
"(" << i+1
1636 <<
") = " << (*(unifiedHistQuanttsForAllBins[j]))[i]
1646 "BaseVectorSequence<V,M>::computeHistCdfstaccKde()",
1647 "unified histogram writing, parallel vectors not supported yet");
1651 for (
unsigned int i = 0; i < unifiedHistQuanttsForAllBins.size(); ++i) {
1652 if (unifiedHistQuanttsForAllBins[i] != NULL)
delete unifiedHistQuanttsForAllBins[i];
1654 for (
unsigned int i = 0; i < unifiedHistCentersForAllBins.size(); ++i) {
1655 if (unifiedHistCentersForAllBins[i] != NULL)
delete unifiedHistCentersForAllBins[i];
1661 if (m_env.subDisplayFile()) {
1662 *m_env.subDisplayFile() <<
"Chain histograms took " << tmpRunTime
1668 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
1672 if ((statisticalOptions.cdfStaccCompute() ) &&
1673 (statisticalOptions.cdfStaccNumEvalPositions() > 0)) {
1675 iRC = gettimeofday(&timevalTmp, NULL);
1676 if (m_env.subDisplayFile()) {
1677 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1678 <<
"\nComputing cdf statistical accuracy"
1682 std::vector<V*> cdfStaccEvalPositions(statisticalOptions.cdfStaccNumEvalPositions(),NULL);
1685 cdfStaccEvalPositions);
1687 std::vector<V*> cdfStaccValues(statisticalOptions.cdfStaccNumEvalPositions(),NULL);
1688 this->subCdfStacc(0,
1689 cdfStaccEvalPositions,
1694 if (m_env.subDisplayFile()) {
1695 *m_env.subDisplayFile() <<
"Chain cdf statistical accuracy took " << tmpRunTime
1704 if ((statisticalOptions.kdeCompute() ) &&
1705 (statisticalOptions.kdeNumEvalPositions() > 0)) {
1707 iRC = gettimeofday(&timevalTmp, NULL);
1708 if (m_env.subDisplayFile()) {
1709 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1710 <<
"\nComputing Kde"
1714 std::string subCoreName_GaussianKdePositions((std::string)(
"_GkdePosits_sub")+m_env.subIdString());
1715 std::string uniCoreName_GaussianKdePositions((std::string)(
"_unifGkdePosits_sub")+m_env.subIdString());
1716 if (m_env.numSubEnvironments() == 1) subCoreName_GaussianKdePositions = uniCoreName_GaussianKdePositions;
1718 std::string subCoreName_GaussianKdeScaleVec ((std::string)(
"_GkdeScalev_sub")+m_env.subIdString());
1719 std::string uniCoreName_GaussianKdeScaleVec ((std::string)(
"_unifGkdeScalev_sub")+m_env.subIdString());
1720 if (m_env.numSubEnvironments() == 1) subCoreName_GaussianKdeScaleVec = uniCoreName_GaussianKdeScaleVec;
1722 std::string subCoreName_GaussianKdeValues ((std::string)(
"_GkdeValues_sub")+m_env.subIdString());
1723 std::string uniCoreName_GaussianKdeValues ((std::string)(
"_unifGkdeValues_sub")+m_env.subIdString());
1724 if (m_env.numSubEnvironments() == 1) subCoreName_GaussianKdeValues = uniCoreName_GaussianKdeValues;
1726 V iqrVec(m_vectorSpace.zeroVector());
1727 this->subInterQuantileRange(0,
1730 V gaussianKdeScaleVec(m_vectorSpace.zeroVector());
1731 this->subScalesForKde(0,
1734 gaussianKdeScaleVec);
1736 std::vector<V*> kdeEvalPositions(statisticalOptions.kdeNumEvalPositions(),NULL);
1741 std::vector<V*> gaussianKdeDensities(statisticalOptions.kdeNumEvalPositions(),NULL);
1742 this->subGaussian1dKde(0,
1743 gaussianKdeScaleVec,
1745 gaussianKdeDensities);
1748 if (m_env.subDisplayFile()) {
1749 *m_env.subDisplayFile() <<
"\nComputed inter quantile ranges for chain '" << m_name <<
"'"
1755 *m_env.subDisplayFile() << line;
1757 sprintf(line,
"%9s%s",
1760 *m_env.subDisplayFile() << line;
1762 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1763 sprintf(line,
"\n%8.8s",
1764 m_vectorSpace.localComponentName(i).c_str() );
1765 *m_env.subDisplayFile() << line;
1767 sprintf(line,
"%2s%11.4e",
1770 *m_env.subDisplayFile() << line;
1772 *m_env.subDisplayFile() << std::endl;
1776 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() > 10)) {
1777 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeHistCdfstaccKde()"
1778 <<
", for chain '" << m_name <<
"'"
1779 <<
", about to write sub kde to ofstream with pointer = " << passedOfs
1783 std::ofstream& ofsvar = *passedOfs;
1784 ofsvar << m_name << subCoreName_GaussianKdePositions <<
" = zeros(" << this->vectorSizeLocal()
1785 <<
"," << kdeEvalPositions.size()
1788 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1789 for (
unsigned int j = 0; j < kdeEvalPositions.size(); ++j) {
1790 ofsvar << m_name << subCoreName_GaussianKdePositions <<
"(" << i+1
1792 <<
") = " << (*(kdeEvalPositions[j]))[i]
1798 ofsvar << m_name << subCoreName_GaussianKdeScaleVec <<
" = zeros(" << this->vectorSizeLocal()
1801 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1802 ofsvar << m_name << subCoreName_GaussianKdeScaleVec <<
"(" << i+1
1803 <<
") = " << gaussianKdeScaleVec[i]
1808 ofsvar << m_name << subCoreName_GaussianKdeValues <<
" = zeros(" << this->vectorSizeLocal()
1809 <<
"," << gaussianKdeDensities.size()
1812 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1813 for (
unsigned int j = 0; j < gaussianKdeDensities.size(); ++j) {
1814 ofsvar << m_name << subCoreName_GaussianKdeValues <<
"(" << i+1
1816 <<
") = " << (*(gaussianKdeDensities[j]))[i]
1823 for (
unsigned int i = 0; i < gaussianKdeDensities.size(); ++i) {
1824 if (gaussianKdeDensities[i] != NULL)
delete gaussianKdeDensities[i];
1826 for (
unsigned int i = 0; i < kdeEvalPositions.size(); ++i) {
1827 if (kdeEvalPositions[i] != NULL)
delete kdeEvalPositions[i];
1830 if ((
int) m_env.numSubEnvironments() > 1) {
1832 V unifiedIqrVec(m_vectorSpace.zeroVector());
1833 this->unifiedInterQuantileRange(0,
1837 V unifiedGaussianKdeScaleVec(m_vectorSpace.zeroVector());
1838 this->unifiedScalesForKde(0,
1841 unifiedGaussianKdeScaleVec);
1844 std::vector<V*> unifiedKdeEvalPositions(statisticalOptions.kdeNumEvalPositions(),NULL);
1846 unifiedStatsMaxPositions,
1847 unifiedKdeEvalPositions);
1849 std::vector<V*> unifiedGaussianKdeDensities(statisticalOptions.kdeNumEvalPositions(),NULL);
1850 this->unifiedGaussian1dKde(0,
1851 unifiedGaussianKdeScaleVec,
1852 unifiedKdeEvalPositions,
1853 unifiedGaussianKdeDensities);
1857 if (m_env.subDisplayFile()) {
1858 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1859 if (m_env.inter0Rank() == 0) {
1860 *m_env.subDisplayFile() <<
"\nComputed unified inter quantile ranges for chain '" << m_name <<
"'"
1866 *m_env.subDisplayFile() << line;
1868 sprintf(line,
"%9s%s",
1871 *m_env.subDisplayFile() << line;
1873 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1874 sprintf(line,
"\n%8.8s",
1875 m_vectorSpace.localComponentName(i).c_str() );
1876 *m_env.subDisplayFile() << line;
1878 sprintf(line,
"%2s%11.4e",
1881 *m_env.subDisplayFile() << line;
1883 *m_env.subDisplayFile() << std::endl;
1889 "BaseVectorSequence<V,M>::computeHistCdfstaccKde()",
1890 "unified iqr writing, parallel vectors not supported yet");
1895 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() > 10)) {
1896 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeHistCdfstaccKde()"
1897 <<
", for chain '" << m_name <<
"'"
1898 <<
", about to write unified kde to ofstream with pointer = " << passedOfs
1902 if (m_vectorSpace.numOfProcsForStorage() == 1) {
1903 if (m_env.inter0Rank() == 0) {
1904 std::ofstream& ofsvar = *passedOfs;
1905 ofsvar << m_name << uniCoreName_GaussianKdePositions <<
" = zeros(" << this->vectorSizeLocal()
1906 <<
"," << unifiedKdeEvalPositions.size()
1909 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() > 10)) {
1910 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeHistCdfstaccKde()"
1911 <<
", for chain '" << m_name <<
"'"
1912 <<
": just wrote '... = zeros(.,.);' line to output file, which has pointer = " << passedOfs
1915 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1916 for (
unsigned int j = 0; j < unifiedKdeEvalPositions.size(); ++j) {
1917 ofsvar << m_name << uniCoreName_GaussianKdePositions <<
"(" << i+1
1919 <<
") = " << (*(unifiedKdeEvalPositions[j]))[i]
1925 ofsvar << m_name << uniCoreName_GaussianKdeScaleVec <<
" = zeros(" << this->vectorSizeLocal()
1928 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1929 ofsvar << m_name << uniCoreName_GaussianKdeScaleVec <<
"(" << i+1
1930 <<
") = " << unifiedGaussianKdeScaleVec[i]
1935 ofsvar << m_name << uniCoreName_GaussianKdeValues <<
" = zeros(" << this->vectorSizeLocal()
1936 <<
"," << unifiedGaussianKdeDensities.size()
1939 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
1940 for (
unsigned int j = 0; j < unifiedGaussianKdeDensities.size(); ++j) {
1941 ofsvar << m_name << uniCoreName_GaussianKdeValues <<
"(" << i+1
1943 <<
") = " << (*(unifiedGaussianKdeDensities[j]))[i]
1953 "BaseVectorSequence<V,M>::computeHistCdfstaccKde()",
1954 "unified Kde writing, parallel vectors not supported yet");
1958 for (
unsigned int i = 0; i < unifiedGaussianKdeDensities.size(); ++i) {
1959 if (unifiedGaussianKdeDensities[i] != NULL)
delete unifiedGaussianKdeDensities[i];
1961 for (
unsigned int i = 0; i < unifiedKdeEvalPositions.size(); ++i) {
1962 if (unifiedKdeEvalPositions[i] != NULL)
delete unifiedKdeEvalPositions[i];
1968 if (m_env.subDisplayFile()) {
1969 *m_env.subDisplayFile() <<
"Chain Kde took " << tmpRunTime
1975 if (m_env.subDisplayFile()) {
1976 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
1977 <<
"\n Finished computing histogram and/or cdf stacc and/or Kde for chain '" << m_name <<
"'"
1978 <<
"\n-----------------------------------------------------"
1986 template<
class V,
class M>
1988 BaseVectorSequence<V,M>::computeCovCorrMatrices(
1989 const SequenceStatisticalOptions& statisticalOptions,
1990 std::ofstream* passedOfs)
1992 if (m_env.subDisplayFile()) {
1993 *m_env.subDisplayFile() <<
"\n"
1994 <<
"\n-----------------------------------------------------"
1995 <<
"\n Computing covariance and correlation matrices for chain '" << m_name <<
"' ..."
1996 <<
"\n-----------------------------------------------------"
2003 M* covarianceMatrix =
new M(m_env,
2004 m_vectorSpace.map(),
2005 m_vectorSpace.dimGlobal());
2006 M* correlationMatrix =
new M(m_env,
2007 m_vectorSpace.map(),
2008 m_vectorSpace.dimGlobal());
2012 this->subSequenceSize(),
2014 *correlationMatrix);
2016 if (m_env.subDisplayFile()) {
2017 if (m_vectorSpace.numOfProcsForStorage() == 1) {
2019 if (m_env.inter0Rank() == 0) {
2020 *m_env.subDisplayFile() <<
"\nBaseVectorSequence<V,M>::computeCovCorrMatrices"
2021 <<
", chain " << m_name
2022 <<
": contents of covariance matrix are\n" << *covarianceMatrix
2025 *m_env.subDisplayFile() <<
"\nBaseVectorSequence<V,M>::computeCovCorrMatrices"
2026 <<
", chain " << m_name
2027 <<
": contents of correlation matrix are\n" << *correlationMatrix
2034 "BaseVectorSequence<V,M>::computeCovCorrMatrices()",
2035 "parallel vectors not supported yet");
2039 delete correlationMatrix;
2040 delete covarianceMatrix;
2042 if (m_env.subDisplayFile()) {
2043 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2044 <<
"\n Finished computing covariance and correlation matrices for chain '" << m_name <<
"'"
2045 <<
"\n-----------------------------------------------------"
2052 #endif // #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2058 #ifdef UQ_CODE_HAS_MONITORS
2059 template<
class V,
class M>
2061 BaseVectorSequence<V,M>::computeMeanEvolution(
2062 const SequenceStatisticalOptions& statisticalOptions,
2063 std::ofstream* passedOfs)
2066 struct timeval timevalTmp;
2067 iRC = gettimeofday(&timevalTmp, NULL);
2068 double tmpRunTime = 0.;
2070 if (m_env.subDisplayFile()) {
2071 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2072 <<
"\nComputing mean evolution"
2076 unsigned int monitorPeriod = statisticalOptions.meanMonitorPeriod();
2077 unsigned int iMin = 0;
2078 unsigned int iMax = (
unsigned int) ( ((
double) this->subSequenceSize())/((double) monitorPeriod) );
2080 if (monitorPeriod == 1) {
2084 if (m_env.subDisplayFile()) {
2085 *m_env.subDisplayFile() <<
" Sub sequence size = " << this->subSequenceSize()
2086 <<
"\n Monitor period = " << monitorPeriod
2087 <<
"\n Number of monitoring positions = " << iMax
2091 this->subMeanMonitorAlloc(iMax);
2092 if (m_env.numSubEnvironments() > 1) {
2093 this->subMeanInter0MonitorAlloc(iMax);
2094 this->unifiedMeanMonitorAlloc(iMax);
2097 for (
unsigned int i = iMin; i < iMax; ++i) {
2098 unsigned int currentMonitoredFinalPosition = (i+1)*monitorPeriod - 1;
2099 V subMeanVec (m_vectorSpace.zeroVector());
2100 V subMeanCltStd(m_vectorSpace.zeroVector());
2101 this->subMeanMonitorRun(currentMonitoredFinalPosition,
2104 this->subMeanMonitorStore(i,
2105 currentMonitoredFinalPosition,
2109 if (m_env.numSubEnvironments() > 1) {
2110 V subMeanInter0Mean (m_vectorSpace.zeroVector());
2111 V subMeanInter0Clt95 (m_vectorSpace.zeroVector());
2112 V subMeanInter0Empirical90(m_vectorSpace.zeroVector());
2113 V subMeanInter0Min (m_vectorSpace.zeroVector());
2114 V subMeanInter0Max (m_vectorSpace.zeroVector());
2115 this->subMeanInter0MonitorRun(currentMonitoredFinalPosition,
2118 subMeanInter0Empirical90,
2121 this->subMeanInter0MonitorStore(i,
2122 currentMonitoredFinalPosition,
2125 subMeanInter0Empirical90,
2129 V unifiedMeanVec (m_vectorSpace.zeroVector());
2130 V unifiedMeanCltStd(m_vectorSpace.zeroVector());
2131 this->unifiedMeanMonitorRun(currentMonitoredFinalPosition,
2134 this->unifiedMeanMonitorStore(i,
2135 currentMonitoredFinalPosition,
2142 this->subMeanMonitorWrite(*passedOfs);
2143 if (m_env.numSubEnvironments() > 1) {
2144 this->subMeanInter0MonitorWrite(*passedOfs);
2145 this->unifiedMeanMonitorWrite(*passedOfs);
2149 this->subMeanMonitorFree();
2150 if (m_env.numSubEnvironments() > 1) {
2151 this->subMeanInter0MonitorFree();
2152 this->unifiedMeanMonitorFree();
2156 if (m_env.subDisplayFile()) {
2157 *m_env.subDisplayFile() <<
"Mean evolution took " << tmpRunTime
2164 #endif //#ifdef UQ_CODE_HAS_MONITORS
2171 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2172 template<
class V,
class M>
2174 BaseVectorSequence<V,M>::computeBMM(
2175 const SequenceStatisticalOptions& statisticalOptions,
2176 const std::vector<unsigned int>& initialPosForStatistics,
2177 std::ofstream* passedOfs)
2180 struct timeval timevalTmp;
2181 iRC = gettimeofday(&timevalTmp, NULL);
2182 double tmpRunTime = 0.;
2184 if (m_env.subDisplayFile()) {
2185 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2186 <<
"\nComputing variance of sample mean through BMM"
2190 if (m_env.subDisplayFile()) {
2191 *m_env.subDisplayFile() <<
"In BaseVectorSequence<V,M>::computeBMM(): lengths for batchs in BMM =";
2192 for (
unsigned int i = 0; i < statisticalOptions.bmmLengths().size(); ++i) {
2193 *m_env.subDisplayFile() <<
" " << statisticalOptions.bmmLengths()[i];
2195 *m_env.subDisplayFile() << std::endl;
2198 TwoDArray<V> _2dArrayOfBMM(initialPosForStatistics.size(),statisticalOptions.bmmLengths().size());
2199 for (
unsigned int i = 0; i < _2dArrayOfBMM.numRows(); ++i) {
2200 for (
unsigned int j = 0; j < _2dArrayOfBMM.numCols(); ++j) {
2201 _2dArrayOfBMM.setLocation(i,j,
new V(m_vectorSpace.zeroVector()) );
2204 V bmmVec(m_vectorSpace.zeroVector());
2205 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2206 unsigned int initialPos = initialPosForStatistics[initialPosId];
2207 for (
unsigned int batchLengthId = 0; batchLengthId < statisticalOptions.bmmLengths().size(); batchLengthId++) {
2208 unsigned int batchLength = statisticalOptions.bmmLengths()[batchLengthId];
2209 this->bmm(initialPos,
2212 _2dArrayOfBMM(initialPosId,batchLengthId) = bmmVec;
2216 if (m_env.subDisplayFile()) {
2217 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2218 *m_env.subDisplayFile() <<
"\nEstimated variance of sample mean, through batch means method, for subchain beginning at position " << initialPosForStatistics[initialPosId]
2219 <<
" (each column corresponds to a batch length)"
2225 *m_env.subDisplayFile() << line;
2226 for (
unsigned int batchLengthId = 0; batchLengthId < statisticalOptions.bmmLengths().size(); batchLengthId++) {
2227 sprintf(line,
"%10s%3d",
2229 statisticalOptions.bmmLengths()[batchLengthId]);
2230 *m_env.subDisplayFile() << line;
2233 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2234 sprintf(line,
"\n%9.9s",
2235 m_vectorSpace.localComponentName(i).c_str() );
2236 *m_env.subDisplayFile() << line;
2237 for (
unsigned int batchLengthId = 0; batchLengthId < statisticalOptions.bmmLengths().size(); batchLengthId++) {
2238 sprintf(line,
"%2s%11.4e",
2240 _2dArrayOfBMM(initialPosId,batchLengthId)[i]);
2241 *m_env.subDisplayFile() << line;
2244 *m_env.subDisplayFile() << std::endl;
2249 if (m_env.subDisplayFile()) {
2250 *m_env.subDisplayFile() <<
"Chain BMM took " << tmpRunTime
2258 template<
class V,
class M>
2260 BaseVectorSequence<V,M>::computeFFT(
2261 const SequenceStatisticalOptions& statisticalOptions,
2262 const std::vector<unsigned int>& initialPosForStatistics,
2263 std::ofstream* passedOfs)
2266 struct timeval timevalTmp;
2267 iRC = gettimeofday(&timevalTmp, NULL);
2268 double tmpRunTime = 0.;
2270 if (m_env.subDisplayFile()) {
2271 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2272 <<
"\nComputing FFT of chain on parameter of id = " << statisticalOptions.fftParamId()
2276 std::vector<std::complex<double> > forwardResult(0,std::complex<double>(0.,0.));
2277 std::vector<std::complex<double> > inverseResult(0,std::complex<double>(0.,0.));
2278 Fft<std::complex<double> > fftObj(m_env);
2279 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2280 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2281 this->fftForward(initialPosition,
2282 statisticalOptions.fftSize(),
2283 statisticalOptions.fftParamId(),
2286 if (statisticalOptions.fftWrite() && passedOfs) {
2287 std::ofstream& ofsvar = *passedOfs;
2288 ofsvar << m_name <<
"_fftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << 1
2289 <<
"," << forwardResult.size()
2292 for (
unsigned int j = 0; j < forwardResult.size(); ++j) {
2293 ofsvar << m_name <<
"_fftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << 1
2295 <<
") = " << forwardResult[j].real()
2296 <<
" + i*" << forwardResult[j].imag()
2302 if (statisticalOptions.fftTestInversion()) {
2303 fftObj.inverse(forwardResult,
2304 statisticalOptions.fftSize(),
2306 if (statisticalOptions.fftWrite() && passedOfs) {
2307 std::ofstream& ofsvar = *passedOfs;
2308 ofsvar << m_name <<
"_iftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << 1
2309 <<
"," << inverseResult.size()
2312 for (
unsigned int j = 0; j < inverseResult.size(); ++j) {
2313 ofsvar << m_name <<
"_iftInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << 1
2315 <<
") = " << inverseResult[j].real()
2316 <<
" + i*" << inverseResult[j].imag()
2325 if (m_env.subDisplayFile()) {
2326 *m_env.subDisplayFile() <<
"Chain FFT took " << tmpRunTime
2334 template<
class V,
class M>
2336 BaseVectorSequence<V,M>::computePSD(
2337 const SequenceStatisticalOptions& statisticalOptions,
2338 const std::vector<unsigned int>& initialPosForStatistics,
2339 std::ofstream* passedOfs)
2342 struct timeval timevalTmp;
2343 iRC = gettimeofday(&timevalTmp, NULL);
2344 double tmpRunTime = 0.;
2346 if (m_env.subDisplayFile()) {
2347 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2348 <<
"\nComputing PSD of chain on parameter of id = " << statisticalOptions.psdParamId()
2352 std::vector<double> psdResult(0,0.);
2353 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2354 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2355 this->psd(initialPosition,
2356 statisticalOptions.psdNumBlocks(),
2357 statisticalOptions.psdHopSizeRatio(),
2358 statisticalOptions.psdParamId(),
2361 if (statisticalOptions.psdWrite() && passedOfs) {
2362 std::ofstream& ofsvar = *passedOfs;
2363 ofsvar << m_name <<
"_psdInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << 1
2364 <<
"," << psdResult.size()
2367 for (
unsigned int j = 0; j < psdResult.size(); ++j) {
2368 ofsvar << m_name <<
"_psdInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << 1
2370 <<
") = " << psdResult[j]
2378 if (m_env.subDisplayFile()) {
2379 *m_env.subDisplayFile() <<
"Chain PSD took " << tmpRunTime
2387 template<
class V,
class M>
2389 BaseVectorSequence<V,M>::computePSDAtZero(
2390 const SequenceStatisticalOptions& statisticalOptions,
2391 const std::vector<unsigned int>& initialPosForStatistics,
2392 std::ofstream* passedOfs)
2395 struct timeval timevalTmp;
2396 iRC = gettimeofday(&timevalTmp, NULL);
2397 double tmpRunTime = 0.;
2399 if (m_env.subDisplayFile()) {
2400 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2401 <<
"\nComputing PSD at frequency zero for all parameters"
2405 TwoDArray<V> _2dArrayOfPSDAtZero(initialPosForStatistics.size(),statisticalOptions.psdAtZeroNumBlocks().size());
2406 for (
unsigned int i = 0; i < _2dArrayOfPSDAtZero.numRows(); ++i) {
2407 for (
unsigned int j = 0; j < _2dArrayOfPSDAtZero.numCols(); ++j) {
2408 _2dArrayOfPSDAtZero.setLocation(i,j,
new V(m_vectorSpace.zeroVector()) );
2411 V psdVec(m_vectorSpace.zeroVector());
2412 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2413 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2414 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2415 unsigned int numBlocks = statisticalOptions.psdAtZeroNumBlocks()[numBlocksId];
2416 this->psdAtZero(initialPosition,
2418 statisticalOptions.psdAtZeroHopSizeRatio(),
2420 _2dArrayOfPSDAtZero(initialPosId,numBlocksId) = psdVec;
2425 if ((statisticalOptions.psdAtZeroDisplay()) && (m_env.subDisplayFile())) {
2426 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2427 unsigned int initialPos = initialPosForStatistics[initialPosId];
2428 *m_env.subDisplayFile() <<
"\nComputed PSD at frequency zero for subchain beginning at position " << initialPos
2429 <<
", so effective data size = " << this->subSequenceSize() - initialPos
2430 <<
" (each column corresponds to a number of blocks)"
2436 *m_env.subDisplayFile() << line;
2437 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2438 sprintf(line,
"%10s%3d",
2440 statisticalOptions.psdAtZeroNumBlocks()[numBlocksId]);
2441 *m_env.subDisplayFile() << line;
2444 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2445 sprintf(line,
"\n%9.9s",
2446 m_vectorSpace.localComponentName(i).c_str() );
2447 *m_env.subDisplayFile() << line;
2448 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2449 sprintf(line,
"%2s%11.4e",
2451 _2dArrayOfPSDAtZero(initialPosId,numBlocksId)[i]);
2452 *m_env.subDisplayFile() << line;
2455 *m_env.subDisplayFile() << std::endl;
2460 if ( (m_env.subDisplayFile())) {
2461 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2462 unsigned int initialPos = initialPosForStatistics[initialPosId];
2463 *m_env.subDisplayFile() <<
"\nEstimated variance of sample mean, through psd, for subchain beginning at position " << initialPos
2464 <<
", so effective data size = " << this->subSequenceSize() - initialPos
2465 <<
" (each column corresponds to a number of blocks)"
2471 *m_env.subDisplayFile() << line;
2472 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2473 sprintf(line,
"%10s%3d",
2475 statisticalOptions.psdAtZeroNumBlocks()[numBlocksId]);
2476 *m_env.subDisplayFile() << line;
2479 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2480 sprintf(line,
"\n%9.9s",
2481 m_vectorSpace.localComponentName(i).c_str() );
2482 *m_env.subDisplayFile() << line;
2483 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2484 sprintf(line,
"%2s%11.4e",
2486 2.*M_PI*_2dArrayOfPSDAtZero(initialPosId,numBlocksId)[i]/(
double) (this->subSequenceSize() - initialPos));
2487 *m_env.subDisplayFile() << line;
2490 *m_env.subDisplayFile() << std::endl;
2495 if (m_env.subDisplayFile()) {
2496 *m_env.subDisplayFile() <<
"Chain PSD at frequency zero took " << tmpRunTime
2502 if (statisticalOptions.psdAtZeroWrite() && passedOfs) {
2503 std::ofstream& ofsvar = *passedOfs;
2504 ofsvar << m_name <<
"_psdAtZeroNumBlocks_sub" << m_env.subIdString() <<
" = zeros(" << 1
2505 <<
"," << statisticalOptions.psdAtZeroNumBlocks().size()
2508 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2509 ofsvar << m_name <<
"_psdAtZeroNumBlocks_sub" << m_env.subIdString() <<
"(" << 1
2510 <<
"," << numBlocksId+1
2511 <<
") = " << statisticalOptions.psdAtZeroNumBlocks()[numBlocksId]
2516 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2517 ofsvar << m_name <<
"_psdAtZeroInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
" = zeros(" << this->vectorSizeLocal()
2518 <<
"," << statisticalOptions.psdAtZeroNumBlocks().size()
2521 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2522 for (
unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2523 ofsvar << m_name <<
"_psdAtZeroInitPos" << initialPosForStatistics[initialPosId] <<
"_sub" << m_env.subIdString() <<
"(" << i+1
2524 <<
"," << numBlocksId+1
2525 <<
") = " << _2dArrayOfPSDAtZero(initialPosId,numBlocksId)[i]
2536 template<
class V,
class M>
2538 BaseVectorSequence<V,M>::computeGeweke(
2539 const SequenceStatisticalOptions& statisticalOptions,
2540 const std::vector<unsigned int>& initialPosForStatistics,
2541 std::ofstream* passedOfs)
2544 struct timeval timevalTmp;
2545 iRC = gettimeofday(&timevalTmp, NULL);
2546 double tmpRunTime = 0.;
2548 if (m_env.subDisplayFile()) {
2549 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2550 <<
"\nComputing Geweke coefficients"
2554 std::vector<V*> vectorOfGeweke(initialPosForStatistics.size(),NULL);
2555 V gewVec(m_vectorSpace.zeroVector());
2556 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2557 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2558 this->geweke(initialPosition,
2559 statisticalOptions.gewekeNaRatio(),
2560 statisticalOptions.gewekeNbRatio(),
2562 vectorOfGeweke[initialPosId] =
new V(gewVec);
2565 if (m_env.subDisplayFile()) {
2566 *m_env.subDisplayFile() <<
"\nComputed Geweke coefficients with 10% and 50% percentages"
2567 <<
" (each column corresponds to a different initial position on the full chain)"
2573 *m_env.subDisplayFile() << line;
2574 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2575 sprintf(line,
"%10s%3d",
2577 initialPosForStatistics[initialPosId]);
2578 *m_env.subDisplayFile() << line;
2581 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2582 sprintf(line,
"\n%9.9s",
2583 m_vectorSpace.localComponentName(i).c_str() );
2584 *m_env.subDisplayFile() << line;
2585 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2586 sprintf(line,
"%2s%11.4e",
2588 (*(vectorOfGeweke[initialPosId]))[i]);
2589 *m_env.subDisplayFile() << line;
2592 *m_env.subDisplayFile() << std::endl;
2596 if (m_env.subDisplayFile()) {
2597 *m_env.subDisplayFile() <<
"Chain Geweke took " << tmpRunTime
2605 template<
class V,
class M>
2607 BaseVectorSequence<V,M>::computeMeanStacc(
2608 const SequenceStatisticalOptions& statisticalOptions,
2609 const std::vector<unsigned int>& initialPosForStatistics,
2610 std::ofstream* passedOfs)
2613 struct timeval timevalTmp;
2614 iRC = gettimeofday(&timevalTmp, NULL);
2615 double tmpRunTime = 0.;
2617 if (m_env.subDisplayFile()) {
2618 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
2619 <<
"\nComputing mean statistical accuracy"
2623 std::vector<V*> vectorOfMeanStacc(initialPosForStatistics.size(),NULL);
2624 V meanStaccVec(m_vectorSpace.zeroVector());
2625 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2626 unsigned int initialPosition = initialPosForStatistics[initialPosId];
2627 this->meanStacc(initialPosition,
2629 vectorOfMeanStacc[initialPosId] =
new V(meanStaccVec);
2632 if (m_env.subDisplayFile()) {
2633 *m_env.subDisplayFile() <<
"\nComputed mean statistical accuracy"
2634 <<
" (each column corresponds to a different initial position on the full chain)"
2640 *m_env.subDisplayFile() << line;
2641 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2642 sprintf(line,
"%10s%3d",
2644 initialPosForStatistics[initialPosId]);
2645 *m_env.subDisplayFile() << line;
2648 for (
unsigned int i = 0; i < this->vectorSizeLocal() ; ++i) {
2649 sprintf(line,
"\n%9.9s",
2650 m_vectorSpace.localComponentName(i).c_str() );
2651 *m_env.subDisplayFile() << line;
2652 for (
unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2653 sprintf(line,
"%2s%11.4e",
2655 (*(vectorOfMeanStacc[initialPosId]))[i]);
2656 *m_env.subDisplayFile() << line;
2659 *m_env.subDisplayFile() << std::endl;
2663 if (m_env.subDisplayFile()) {
2664 *m_env.subDisplayFile() <<
"Chain mean statistical accuracy took " << tmpRunTime
2671 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2677 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
2682 unsigned int subNumSamples,
2694 "ComputeCovCorrMatricesBetweenVectorSequences()",
2695 "parallel vectors not supported yet");
2700 UQ_FATAL_TEST_MACRO((numRowsLocal != pqCovMatrix.numRowsLocal()) || (numCols != pqCovMatrix.numCols()),
2702 "ComputeCovCorrMatricesBetweenVectorSequences()",
2703 "inconsistent dimensions for covariance matrix");
2705 UQ_FATAL_TEST_MACRO((numRowsLocal != pqCorrMatrix.numRowsLocal()) || (numCols != pqCorrMatrix.numCols()),
2707 "ComputeCorrelationBetweenVectorSequences()",
2708 "inconsistent dimensions for correlation matrix");
2712 "ComputeCovCorrMatricesBetweenVectorSequences()",
2713 "subNumSamples is too large");
2727 for (
unsigned i = 0; i < numRowsLocal; ++i) {
2728 for (
unsigned j = 0; j < numCols; ++j) {
2729 pqCovMatrix(i,j) = 0.;
2732 for (
unsigned k = 0; k < subNumSamples; ++k) {
2735 tmpP -= unifiedMeanP;
2738 tmpQ -= unifiedMeanQ;
2740 for (
unsigned i = 0; i < numRowsLocal; ++i) {
2741 for (
unsigned j = 0; j < numCols; ++j) {
2742 pqCovMatrix(i,j) += tmpP[i]*tmpQ[j];
2752 unifiedSampleVarianceP);
2758 unifiedSampleVarianceQ);
2761 if (useOnlyInter0Comm) {
2763 unsigned int unifiedNumSamples = 0;
2765 "ComputeCovCorrMatricesBetweenVectorSequences()",
2766 "failed MPI.Allreduce() for subNumSamples");
2768 for (
unsigned i = 0; i < numRowsLocal; ++i) {
2769 for (
unsigned j = 0; j < numCols; ++j) {
2772 "ComputeCovCorrMatricesBetweenVectorSequences()",
2773 "failed MPI.Allreduce() for a matrix position");
2774 pqCovMatrix(i,j) = aux/((double) (unifiedNumSamples-1));
2778 for (
unsigned i = 0; i < numRowsLocal; ++i) {
2779 for (
unsigned j = 0; j < numCols; ++j) {
2780 pqCorrMatrix(i,j) = pqCovMatrix(i,j)/std::sqrt(unifiedSampleVarianceP[i])/std::sqrt(unifiedSampleVarianceQ[j]);
2781 if (((pqCorrMatrix(i,j) + 1.) < -1.e-8) ||
2782 ((pqCorrMatrix(i,j) - 1.) > 1.e-8)) {
2784 std::cerr <<
"In ComputeCovCorrMatricesBetweenVectorSequences()"
2788 <<
", pqCorrMatrix(i,j)+1 = " << pqCorrMatrix(i,j)+1.
2789 <<
", pqCorrMatrix(i,j)-1 = " << pqCorrMatrix(i,j)-1.
2795 ((pqCorrMatrix(i,j) - 1.) > 1.e-8),
2797 "ComputeCovCorrMatricesBetweenVectorSequences()",
2798 "computed correlation is out of range");
2809 "ComputeCovCorrMatricesBetweenVectorSequences()",
2810 "parallel vectors not supported yet (2)");
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
unsigned int vectorSizeGlobal() const
Global dimension (size) of the vector space.
double MiscGetEllapsedSeconds(struct timeval *timeval0)
void deleteStoredVectors()
Deletes all the stored vectors.
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...
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization.
Class for matrix operations using GSL library.
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...
const V & unifiedMinPlain() const
Finds the minimum value of the unified sequence.
const V & subMaxPlain() const
Finds the maximum value of the sub-sequence.
const V & unifiedMeanPlain() const
Finds the mean value of the unified sequence.
int inter0Rank() const
Returns the process inter0 rank.
int worldRank() const
Returns the process world rank.
const V & subSampleVariancePlain() const
Finds the variance of a sample of the sub-sequence.
unsigned int dimLocal() const
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...
#define RawValue_MPI_DOUBLE
#define RawValue_MPI_UNSIGNED
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
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)
unsigned int vectorSizeLocal() const
Local dimension (size) of the vector space.
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.
Class representing a subset of a vector space shaped like a hypercube.
unsigned int numOfProcsForStorage() const
Returns total number of processes.
const T & subMaxPlain() const
Finds the maximum value of the sub-sequence of scalars.
virtual ~BaseVectorSequence()
Destructor.
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
const V & unifiedMaxPlain() const
Finds the maximum value of the unified sequence.
void Allreduce(void *sendbuf, void *recvbuf, int count, RawType_MPI_Datatype datatype, RawType_MPI_Op op, const char *whereMsg, const char *whatMsg) const
Combines values from all processes and distributes the result back to all processes.
const VectorSpace< V, M > & vectorSpace() const
Vector space; access to protected attribute VectorSpace<V,M>& m_vectorSpace.
const V & subMinPlain() const
Finds the minimum value of the sub-sequence.
void copy(const BaseVectorSequence< V, M > &src)
Copies vector sequence src to this.
const V & subMeanPlain() const
Finds the mean value of the sub-sequence.
BaseVectorSequence(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
void append(const BaseVectorSequence< V, M > &src, unsigned int initialPos, unsigned int numPos)
Appends the vector src to this vector.
const V & unifiedMedianPlain() const
Finds the median value of the unified sequence.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
double subPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &subPositionsOfMaximum)
Finds the positions where the maximum element occurs in the sub-sequence.
void MiscComputePositionsBetweenMinMax(V minValues, V maxValues, std::vector< V * > &positions)
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...
Class for a Fast Fourier Transform (FFT) algorithm.
const V & subMedianPlain() const
Finds the median value of the sub-sequence.
const std::string & name() const
Access to protected attribute m_name: name of the sequence of vectors.
Base class for handling vector and array samples (sequence of vectors or arrays). ...
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
A class representing a vector space.
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...
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 MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
const V & unifiedSampleVariancePlain() const
Finds the variance of a sample of the unified sequence.
void computeFilterParams(std::ofstream *passedOfs, unsigned int &initialPos, unsigned int &spacing)
Computes the filtering parameters spacing for the sequence of vectors.
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 ...
double unifiedPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &unifiedPositionsOfMaximum)
Finds the positions where the maximum element occurs in the unified sequence.
unsigned int unifiedSequenceSize() const
Calculates the size of the unified sequence of vectors.
const VectorSpace< V, M > & m_vectorSpace
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
void clear()
Reset the values and the size of the sequence of vectors.
unsigned int dimGlobal() const