Class for handling scalar samples. More...

#include <ScalarSequence.h>

Collaboration diagram for QUESO::ScalarSequence< T >:

Public Types
Class typedefs
typedef std::vector< T >::iterator	seqScalarPositionIteratorTypedef

typedef std::vector< T > ::const_iterator	seqScalarPositionConstIteratorTypedef

Public Member Functions
Constructor/Destructor methods
	ScalarSequence (const BaseEnvironment &env, unsigned int subSequenceSize, const std::string &name)
	Default constructor. More...

	~ScalarSequence ()
	Destructor. More...

Set methods
ScalarSequence< T > &	operator= (const ScalarSequence< T > &rhs)
	Assignment operator; it copies `rhs` to `this`. More...

Access methods
const T &	operator[] (unsigned int posId) const
	Access position `posId` of the sequence of scalars (const). More...

T &	operator[] (unsigned int posId)
	Access position `posId` of the sequence of scalars (non-const). More...

Sequence methods
const BaseEnvironment &	env () const
	Access to QUESO environment. More...

const std::string &	name () const
	Access to the name of the sequence of scalars. More...

void	setName (const std::string &newName)
	Sets a new name to the sequence of scalars. More...

void	clear ()
	Clears the sequence of scalars. More...

unsigned int	subSequenceSize () const
	Size of the sub-sequence of scalars. More...

unsigned int	unifiedSequenceSize (bool useOnlyInter0Comm) const
	Size of the unified sequence of scalars. More...

void	resizeSequence (unsigned int newSequenceSize)
	Resizes the size of the sequence of scalars. More...

void	resetValues (unsigned int initialPos, unsigned int)
	Sets `numPos` values of the sequence to zero, starting at position `initialPos`. More...

void	erasePositions (unsigned int initialPos, unsigned int numPos)
	Erases `numPos` values of the sequence, starting at position `initialPos`. More...

void	getUnifiedContentsAtProc0Only (bool useOnlyInter0Comm, std::vector< T > &outputVec) const
	Gets the unified contents of processor of rank equals to 0. More...

const T &	subMinPlain () const
	Finds the minimum value of the sub-sequence of scalars. More...

const T &	unifiedMinPlain (bool useOnlyInter0Comm) const
	Finds the minimum value of the unified sequence of scalars. More...

const T &	subMaxPlain () const
	Finds the maximum value of the sub-sequence of scalars. More...

const T &	unifiedMaxPlain (bool useOnlyInter0Comm) const
	Finds the maximum value of the unified sequence of scalars. More...

const T &	subMeanPlain () const
	Finds the mean value of the sub-sequence of scalars. More...

const T &	unifiedMeanPlain (bool useOnlyInter0Comm) const
	Finds the mean value of the unified sequence of scalars. More...

const T &	subMedianPlain () const
	Finds the median value of the sub-sequence of scalars. More...

const T &	unifiedMedianPlain (bool useOnlyInter0Comm) const
	Finds the median value of the unified sequence of scalars. More...

const T &	subSampleVariancePlain () const
	Finds the variance of a sample of the sub-sequence of scalars. More...

const T &	unifiedSampleVariancePlain (bool useOnlyInter0Comm) const
	Finds the variance of a sample of the unified sequence of scalars. More...

void	deleteStoredScalars ()
	Deletes all stored scalars. More...

void	setGaussian (const T &mean, const T &stdDev)
	Sets the values of the sequence as a Gaussian distribution of mean given by `meanVec` and standard deviation by `stdDevVec`. More...

void	setUniform (const T &a, const T &b)
	Sets the values of the sequence as a uniform distribution between the values given by vectors `aVec` and `bVec`. More...

void	subUniformlySampledCdf (unsigned int numIntervals, T &minDomainValue, T &maxDomainValue, std::vector< T > &cdfValues) const
	Uniformly samples from the CDF from the sub-sequence. More...

void	unifiedUniformlySampledCdf (bool useOnlyInter0Comm, unsigned int numIntervals, T &unifiedMinDomainValue, T &unifiedMaxDomainValue, std::vector< T > &unifiedCdfValues) const
	Uniformly samples from the CDF from the unified sequence. More...

void	subBasicCdf (unsigned int numIntervals, UniformOneDGrid< T > *&gridValues, std::vector< T > &cdfValues) const
	Finds the Cumulative Distribution Function (CDF) of the sub-sequence of scalars. More...

void	subWeightCdf (unsigned int numIntervals, std::vector< T > &gridValues, std::vector< T > &cdfValues) const
	Finds the Weighted Cumulative Distribution Function (CDF) of the sub-sequence of scalars. More...

void	subWeightCdf (unsigned int numIntervals, UniformOneDGrid< T > *&gridValues, std::vector< T > &cdfValues) const
	Finds the Weighted Cumulative Distribution Function (CDF) of the sub-sequence of scalars. More...

T	subMeanExtra (unsigned int initialPos, unsigned int numPos) const
	Finds the mean value of the sub-sequence, considering `numPos` positions starting at position `initialPos`. More...

T	unifiedMeanExtra (bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos) const
	Finds the mean value of the unified sequence of `numPos` positions starting at position `initialPos`. More...

T	subMedianExtra (unsigned int initialPos, unsigned int numPos) const
	Finds the median value of the sub-sequence, considering `numPos` positions starting at position `initialPos`. More...

T	unifiedMedianExtra (bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos) const
	Finds the median value of the unified sequence, considering `numPos` positions starting at position `initialPos`. More...

T	subSampleVarianceExtra (unsigned int initialPos, unsigned int numPos, const T &meanValue) const
	Finds the sample variance of the sub-sequence, considering `numPos` positions starting at position `initialPos` and of mean `meanVec`. More...

T	unifiedSampleVarianceExtra (bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos, const T &unifiedMeanValue) const
	Finds the sample variance of the unified sequence, considering `numPos` positions starting at position `initialPos` and of mean `meanVec`. More...

T	subSampleStd (unsigned int initialPos, unsigned int numPos, const T &meanValue) const
	Finds the sample standard deviation of the unified sequence, considering `numPos` positions starting at position `initialPos` and of mean `meanValue`. More...

T	unifiedSampleStd (bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos, const T &unifiedMeanValue) const
	Finds the sample standard deviation of the unified sequence, considering `localnumPos` positions starting at position `initialPos` and of mean `unifiedMeanValue`. More...

T	subPopulationVariance (unsigned int initialPos, unsigned int numPos, const T &meanValue) const
	Finds the population variance of the sub-sequence, considering `numPos` positions starting at position `initialPos` and of mean `meanValue`. More...

T	unifiedPopulationVariance (bool useOnlyInter0Comm, unsigned int initialPos, unsigned int numPos, const T &unifiedMeanValue) const
	Finds the population variance of the unified sequence, considering `numPos` positions starting at position `initialPos` and of mean `meanValue`. More...

T	autoCovariance (unsigned int initialPos, unsigned int numPos, const T &meanValue, unsigned int lag) const
	Calculates the autocovariance. More...

T	autoCorrViaDef (unsigned int initialPos, unsigned int numPos, unsigned int lag) const
	Calculates the autocorrelation via definition. More...

void	autoCorrViaFft (unsigned int initialPos, unsigned int numPos, unsigned int maxLag, std::vector< T > &autoCorrs) const
	Calculates the autocorrelation via Fast Fourier transforms (FFT). More...

void	autoCorrViaFft (unsigned int initialPos, unsigned int numPos, unsigned int numSum, T &autoCorrsSum) const
	Calculates the sum of autocorrelation via Fast Fourier transforms (FFT). More...

void	subMinMaxExtra (unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
	Finds the minimum and the maximum values of the sub-sequence, considering `numPos` positions starting at position `initialPos`. More...

void	unifiedMinMaxExtra (bool useOnlyInter0Comm, unsigned int initialPos, unsigned int numPos, T &unifiedMinValue, T &unifiedMaxValue) const
	Finds the minimum and the maximum values of the unified sequence, considering `numPos` positions starting at position `initialPos`. More...

void	subHistogram (unsigned int initialPos, const T &minHorizontalValue, const T &maxHorizontalValue, std::vector< T > &centers, std::vector< unsigned int > &bins) const
	Calculates the histogram of the sub-sequence. More...

void	unifiedHistogram (bool useOnlyInter0Comm, unsigned int initialPos, const T &unifiedMinHorizontalValue, const T &unifiedMaxHorizontalValue, std::vector< T > &unifiedCenters, std::vector< unsigned int > &unifiedBins) const
	Calculates the histogram of the unified sequence. More...

void	subBasicHistogram (unsigned int initialPos, const T &minHorizontalValue, const T &maxHorizontalValue, UniformOneDGrid< T > *&gridValues, std::vector< unsigned int > &bins) const
	Calculates the histogram of the sub-sequence. More...

void	subWeightHistogram (unsigned int initialPos, const T &minHorizontalValue, const T &maxHorizontalValue, UniformOneDGrid< T > *&gridValues, std::vector< unsigned int > &bins) const
	Calculates the weighted histogram of the sub-sequence. More...

void	subWeightHistogram (unsigned int initialPos, const T &minHorizontalValue, const T &maxHorizontalValue, std::vector< T > &gridValues, std::vector< unsigned int > &bins) const
	Calculates the weighted histogram of the sub-sequence. More...

void	subSort (unsigned int initialPos, ScalarSequence< T > &sortedSequence) const
	Sorts the sub-sequence of scalars. More...

void	unifiedSort (bool useOnlyInter0Comm, unsigned int initialPos, ScalarSequence< T > &unifiedSortedSequence) const
	Sorts the unified sequence of scalars. More...

T	subInterQuantileRange (unsigned int initialPos) const
	Returns the interquartile range of the values in the sub-sequence. More...

T	unifiedInterQuantileRange (bool useOnlyInter0Comm, unsigned int initialPos) const
	Returns the interquartile range of the values in the unified sequence. More...

T	subScaleForKde (unsigned int initialPos, const T &iqrValue, unsigned int kdeDimension) const
	Selects the scales (output value) for the kernel density estimation, considering only the sub-sequence. More...

T	unifiedScaleForKde (bool useOnlyInter0Comm, unsigned int initialPos, const T &unifiedIqrValue, unsigned int kdeDimension) const
	Selects the scales (bandwidth) for the kernel density estimation, considering the unified sequence. More...

void	subGaussian1dKde (unsigned int initialPos, double scaleValue, const std::vector< T > &evaluationPositions, std::vector< double > &densityValues) const
	Gaussian kernel for the KDE estimate of the sub-sequence. More...

void	unifiedGaussian1dKde (bool useOnlyInter0Comm, unsigned int initialPos, double unifiedScaleValue, const std::vector< T > &unifiedEvaluationPositions, std::vector< double > &unifiedDensityValues) const
	Gaussian kernel for the KDE estimate of the unified sequence. More...

void	filter (unsigned int initialPos, unsigned int spacing)
	Filters positions in the sequence of vectors. More...

T	brooksGelmanConvMeasure (bool useOnlyInter0Comm, unsigned int initialPos, unsigned int spacing) const
	Estimates convergence rate using Brooks & Gelman method. More...

void	append (const ScalarSequence< T > &src, unsigned int srcInitialPos, unsigned int srcNumPos)
	Appends the scalar sequence `src` to `this` sequence. More...

T	subPositionsOfMaximum (const ScalarSequence< T > &subCorrespondingScalarValues, ScalarSequence< T > &subPositionsOfMaximum)
	Finds the positions where the maximum element occurs in the sub-sequence. More...

T	unifiedPositionsOfMaximum (const ScalarSequence< T > &subCorrespondingScalarValues, ScalarSequence< T > &unifiedPositionsOfMaximum)
	Finds the positions where the maximum element occurs in the unified sequence. More...

void	subWriteContents (unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
	Writes the sub-sequence to a file. More...

void	subWriteContents (unsigned int initialPos, unsigned int numPos, std::ofstream &ofs, const std::string &fileType) const
	Writes the sub-sequence to a file. More...

void	unifiedWriteContents (const std::string &fileName, const std::string &fileType) const
	Writes the unified sequence to a file. More...

void	unifiedReadContents (const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
	Reads the unified sequence from a file. More...

Private Member Functions
void	copy (const ScalarSequence< T > &src)
	Copies the scalar sequence `src` to `this`. More...

void	extractScalarSeq (unsigned int initialPos, unsigned int spacing, unsigned int numPos, ScalarSequence< T > &scalarSeq) const
	Extracts a sequence of scalars. More...

void	extractRawData (unsigned int initialPos, unsigned int spacing, unsigned int numPos, std::vector< double > &rawData) const
	Extracts the raw data. More...

std::vector< T > &	rawData ()
	The sequence of scalars. Access to private attribute `m_seq`. More...

void	subSort ()
	Sorts the sequence of scalars in the private attribute `m_seq`. More...

void	parallelMerge (std::vector< T > &sortedBuffer, const std::vector< T > &leafData, unsigned int treeLevel) const
	Sorts/merges data in parallel using MPI. More...

Private Attributes
const BaseEnvironment &	m_env

std::string	m_name

std::vector< T >	m_seq

T *	m_subMinPlain

T *	m_unifiedMinPlain

T *	m_subMaxPlain

T *	m_unifiedMaxPlain

T *	m_subMeanPlain

T *	m_unifiedMeanPlain

T *	m_subMedianPlain

T *	m_unifiedMedianPlain

T *	m_subSampleVariancePlain

T *	m_unifiedSampleVariancePlain

Detailed Description

template<class T = double>
class QUESO::ScalarSequence< T >

Class for handling scalar samples.

This class handles scalar samples generated by an algorithm, as well as operations that can be carried over them, e.g., calculation of means, correlation and covariance matrices.

Definition at line 54 of file ScalarSequence.h.

Member Typedef Documentation

template<class T = double>

typedef std::vector<T>::const_iterator QUESO::ScalarSequence< T >::seqScalarPositionConstIteratorTypedef

Definition at line 60 of file ScalarSequence.h.

template<class T = double>

typedef std::vector<T>::iterator QUESO::ScalarSequence< T >::seqScalarPositionIteratorTypedef

Definition at line 59 of file ScalarSequence.h.

Constructor & Destructor Documentation

template<class T >

QUESO::ScalarSequence< T >::ScalarSequence	(	const BaseEnvironment &	env,
		unsigned int	subSequenceSize,
		const std::string &	name
	)

Default constructor.

Definition at line 32 of file ScalarSequence.C.

   :
   m_env                       (env),
   m_name                      (name),
   m_seq                       (subSequenceSize,0.),
   m_subMinPlain               (NULL),
   m_unifiedMinPlain           (NULL),
   m_subMaxPlain               (NULL),
   m_unifiedMaxPlain           (NULL),
   m_subMeanPlain              (NULL),
   m_unifiedMeanPlain          (NULL),
   m_subMedianPlain            (NULL),
   m_unifiedMedianPlain        (NULL),
   m_subSampleVariancePlain    (NULL),
   m_unifiedSampleVariancePlain(NULL)
 {
 }

template<class T >

QUESO::ScalarSequence< T >::~ScalarSequence ( )

Destructor.

Definition at line 54 of file ScalarSequence.C.

 {
   deleteStoredScalars();
 }

Member Function Documentation

template<class T>

void QUESO::ScalarSequence< T >::append	(	const ScalarSequence< T > &	src,
		unsigned int	srcInitialPos,
		unsigned int	srcNumPos
	)

Appends the scalar sequence src to this sequence.

This routine deletes all stored computed scalars.

Definition at line 2497 of file ScalarSequence.C.

References QUESO::ScalarSequence< T >::m_seq, queso_require_greater_equal_msg, and QUESO::ScalarSequence< T >::subSequenceSize().

Referenced by QUESO::MLSampling< P_V, P_M >::generateBalLinkedChains_all(), and QUESO::MLSampling< P_V, P_M >::generateUnbLinkedChains_all().

 {
   queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+1), "srcInitialPos is too big");
 
   queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+srcNumPos), "srcNumPos is too big");
 
   deleteStoredScalars();
   unsigned int currentSize = this->subSequenceSize();
   m_seq.resize(currentSize+srcNumPos,0.);
   for (unsigned int i = 0; i < srcNumPos; ++i) {
     m_seq[currentSize+i] = src.m_seq[srcInitialPos+i];
   }
 
   return;
 }

template<class T >

T QUESO::ScalarSequence< T >::autoCorrViaDef	(	unsigned int	initialPos,
		unsigned int	numPos,
		unsigned int	lag
	)		const

Calculates the autocorrelation via definition.

Autocorrelation is the cross-correlation of a variable with itself; it describes the correlation between values of the process at different times, as a function of the two times. It is calculated over a sequence of vectors with initial position initialPos, considering numPos positions, a lag of lag, with mean given by meanValue. Output: the calculated autocorrelations of the sequence of vectors.

Definition at line 1332 of file ScalarSequence.C.

References queso_require_msg.

Referenced by QUESO::SequenceOfVectors< V, M >::autoCorrViaDef().

 {
   bool bRC = ((initialPos          <  this->subSequenceSize()) &&
               (0                   <  numPos                 ) &&
               ((initialPos+numPos) <= this->subSequenceSize()) &&
               (lag                 <  numPos                 )); // lag should not be too large
   queso_require_msg(bRC, "invalid input data");
 
   T meanValue = this->subMeanExtra(initialPos,
                                    numPos);
 
   T covValueZero = this->autoCovariance(initialPos,
                                         numPos,
                                         meanValue,
                                         0); // lag
 
   T corrValue = this->autoCovariance(initialPos,
                                      numPos,
                                      meanValue,
                                      lag);
 
   return corrValue/covValueZero;
 }

template<class T>

void QUESO::ScalarSequence< T >::autoCorrViaFft	(	unsigned int	initialPos,
		unsigned int	numPos,
		unsigned int	maxLag,
		std::vector< T > &	autoCorrs
	)		const

Calculates the autocorrelation via Fast Fourier transforms (FFT).

Definition at line 1361 of file ScalarSequence.C.

References QUESO::Fft< T >::forward(), and QUESO::Fft< T >::inverse().

Referenced by QUESO::ArrayOfSequences< V, M >::autoCorrViaFft(), and QUESO::SequenceOfVectors< V, M >::autoCorrViaFft().

 {
   unsigned int fftSize = 0;
   {
     double tmp = log((double) maxLag)/log(2.);
     double fractionalPart = tmp - ((double) ((unsigned int) tmp));
     if (fractionalPart > 0.) tmp += (1. - fractionalPart);
     unsigned int fftSize1 = (unsigned int) std::pow(2.,tmp+1.); // Yes, tmp+1
     fftSize1 = fftSize1; // To remove warning
 
     tmp = log((double) numPos)/log(2.);
     fractionalPart = tmp - ((double) ((unsigned int) tmp));
     if (fractionalPart > 0.) tmp += (1. - fractionalPart);
     unsigned int fftSize2 = (unsigned int) std::pow(2.,tmp+1);
 
     fftSize = fftSize2;
   }
 
   std::vector<double> rawDataVec(numPos,0.);
   std::vector<std::complex<double> > resultData(0,std::complex<double>(0.,0.));
   Fft<T> fftObj(m_env);
 
   // Forward FFT
   this->extractRawData(initialPos,
                        1, // spacing
                        numPos,
                        rawDataVec);
   T meanValue = this->subMeanExtra(initialPos,
                                    numPos);
   for (unsigned int j = 0; j < numPos; ++j) {
     rawDataVec[j] -= meanValue; // IMPORTANT
   }
 
   rawDataVec.resize(fftSize,0.);
 
   //if (m_env.subDisplayFile()) {
   //  *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
   //                          << ": about to call fftObj.forward()"
   //                          << " with rawDataVec.size() = " << rawDataVec.size()
   //                          << ", fftSize = "            << fftSize
   //                          << ", resultData.size() = "  << resultData.size()
   //                          << std::endl;
   //}
   fftObj.forward(rawDataVec,fftSize,resultData);
 
   // Inverse FFT
   for (unsigned int j = 0; j < fftSize; ++j) {
     rawDataVec[j] = std::norm(resultData[j]);
   }
   //if (m_env.subDisplayFile()) {
   //  *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
   //                          << ": about to call fftObj.inverse()"
   //                          << " with rawDataVec.size() = " << rawDataVec.size()
   //                          << ", fftSize = "            << fftSize
   //                          << ", resultData.size() = "  << resultData.size()
   //                          << std::endl;
   //}
   fftObj.inverse(rawDataVec,fftSize,resultData);
   //if (m_env.subDisplayFile()) {
   //  *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
   //                          << ": returned succesfully from fftObj.inverse()"
   //                          << std::endl;
   //}
 
   // Prepare return data
   autoCorrs.resize(maxLag+1,0.); // Yes, +1
   for (unsigned int j = 0; j < autoCorrs.size(); ++j) {
     double ratio = ((double) j)/((double) (numPos-1));
     autoCorrs[j] = ( resultData[j].real()/resultData[0].real() )*(1.-ratio);
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::autoCorrViaFft	(	unsigned int	initialPos,
		unsigned int	numPos,
		unsigned int	numSum,
		T &	autoCorrsSum
	)		const

Calculates the sum of autocorrelation via Fast Fourier transforms (FFT).

Definition at line 1441 of file ScalarSequence.C.

References QUESO::Fft< T >::forward(), and QUESO::Fft< T >::inverse().

 {
   //if (m_env.subDisplayFile()) {
   //  *m_env.subDisplayFile() << "Entering ScalarSequence<T>::autoCorrViaFft(), for sum"
   //                          << ": initialPos = " << initialPos
   //                          << ", numPos = "     << numPos
   //                          << std::endl;
   //}
 
   double tmp = log((double) numPos)/log(2.);
   double fractionalPart = tmp - ((double) ((unsigned int) tmp));
   if (fractionalPart > 0.) tmp += (1. - fractionalPart);
   unsigned int fftSize = (unsigned int) std::pow(2.,tmp+1);
 
   std::vector<double> rawDataVec(numPos,0.);
   std::vector<std::complex<double> > resultData(0,std::complex<double>(0.,0.));
   Fft<T> fftObj(m_env);
 
   // Forward FFT
   this->extractRawData(initialPos,
                        1, // spacing
                        numPos,
                        rawDataVec);
   T meanValue = this->subMeanExtra(initialPos,
                                    numPos);
   for (unsigned int j = 0; j < numPos; ++j) {
     rawDataVec[j] -= meanValue; // IMPORTANT
   }
   rawDataVec.resize(fftSize,0.);
 
   //if (m_env.subDisplayFile()) {
   //  *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft(), for sum"
   //                          << ": about to call fftObj.forward()"
   //                          << " with rawDataVec.size() = " << rawDataVec.size()
   //                          << ", fftSize = "            << fftSize
   //                          << ", resultData.size() = "  << resultData.size()
   //                          << std::endl;
   //}
   fftObj.forward(rawDataVec,fftSize,resultData);
 
   // Inverse FFT
   for (unsigned int j = 0; j < fftSize; ++j) {
     rawDataVec[j] = std::norm(resultData[j]);
   }
   fftObj.inverse(rawDataVec,fftSize,resultData);
 
   //if (m_env.subDisplayFile()) {
   //  *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft(), for sum"
   //                          << ": computed auto covariance for lag 0 = " << resultData[0].real()/((double) (numPos))
   //                          << ", computed resultData[0].imag() = "      << resultData[0].imag()
   //                          << std::endl;
   //}
 
   // Prepare return data
   autoCorrsSum = 0.;
   for (unsigned int j = 0; j < numSum; ++j) { // Yes, begin at lag '0'
     double ratio = ((double) j)/((double) (numPos-1));
     autoCorrsSum += ( resultData[j].real()/resultData[0].real() )*(1.-ratio);
   }
 
   return;
 }

template<class T>

T QUESO::ScalarSequence< T >::autoCovariance	(	unsigned int	initialPos,
		unsigned int	numPos,
		const T &	meanValue,
		unsigned int	lag
	)		const

Calculates the autocovariance.

The autocovariance is the covariance of a variable with itself at some other time. It is calculated over a sequence of scalars with initial position initialPos, considering numPos positions, a lag of lag, with mean given by meanValue. Output: the calculated autocovariances of the sequence of scalars.

Definition at line 1302 of file ScalarSequence.C.

References queso_require_msg.

Referenced by QUESO::SequenceOfVectors< V, M >::autoCovariance().

 {
   bool bRC = ((initialPos          <  this->subSequenceSize()) &&
               (0                   <  numPos                 ) &&
               ((initialPos+numPos) <= this->subSequenceSize()) &&
               (lag                 <  numPos                 )); // lag should not be too large
   queso_require_msg(bRC, "invalid input data");
 
   unsigned int loopSize      = numPos - lag;
   unsigned int finalPosPlus1 = initialPos + loopSize;
   T diff1;
   T diff2;
   T covValue = 0.;
   for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
     diff1 = m_seq[j    ] - meanValue;
     diff2 = m_seq[j+lag] - meanValue;
     covValue += diff1*diff2;
   }
 
   covValue /= (T) loopSize;
 
   return covValue;
 }

template<class T >

T QUESO::ScalarSequence< T >::brooksGelmanConvMeasure	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		unsigned int	spacing
	)		const

Estimates convergence rate using Brooks & Gelman method.

TODO: implement me!

Definition at line 2468 of file ScalarSequence.C.

References queso_error_msg, and queso_not_implemented.

 {
   double resultValue = 0.;
 
   // FIX ME: put logic if (numSubEnvs == 1) ...
 
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       queso_not_implemented();
     }
     else {
       // Node not in the 'inter0' communicator
       // Do nothing
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return resultValue;
 }

template<class T >

void QUESO::ScalarSequence< T >::clear ( )

Clears the sequence of scalars.

Resets its values and then its size to zero.

Definition at line 123 of file ScalarSequence.C.

Referenced by QUESO::MLSampling< P_V, P_M >::generateSequence_Step02_inter0().

 {
   unsigned int numPos = this->subSequenceSize();
   if (numPos) {
     this->resetValues(0,numPos);
     this->resizeSequence(0);
   }
 
  return;
 }

template<class T>

void QUESO::ScalarSequence< T >::copy ( const ScalarSequence< T > & src )

private

Copies the scalar sequence src to this.

This routine deletes all stored computed scalars.

Definition at line 3102 of file ScalarSequence.C.

References QUESO::ScalarSequence< T >::m_name, QUESO::ScalarSequence< T >::m_seq, and QUESO::ScalarSequence< T >::subSequenceSize().

 {
   m_name = src.m_name;
   m_seq.clear();
   m_seq.resize(src.subSequenceSize(),0.);
   for (unsigned int i = 0; i < m_seq.size(); ++i) {
     m_seq[i] = src.m_seq[i];
   }
   deleteStoredScalars();
 
   return;
 }

template<class T >

void QUESO::ScalarSequence< T >::deleteStoredScalars ( )

Deletes all stored scalars.

Definition at line 465 of file ScalarSequence.C.

 {
   if (m_subMinPlain) {
     delete m_subMinPlain;
     m_subMinPlain                = NULL;
   }
   if (m_unifiedMinPlain) {
     delete m_unifiedMinPlain;
     m_unifiedMinPlain            = NULL;
   }
   if (m_subMaxPlain) {
     delete m_subMaxPlain;
     m_subMaxPlain                = NULL;
   }
   if (m_unifiedMaxPlain) {
     delete m_unifiedMaxPlain;
     m_unifiedMaxPlain            = NULL;
   }
   if (m_subMeanPlain) {
     delete m_subMeanPlain;
     m_subMeanPlain               = NULL;
   }
   if (m_unifiedMeanPlain) {
     delete m_unifiedMeanPlain;
     m_unifiedMeanPlain           = NULL;
   }
   if (m_subMedianPlain) {
     delete m_subMedianPlain;
     m_subMedianPlain             = NULL;
   }
   if (m_unifiedMedianPlain) {
     delete m_unifiedMedianPlain;
     m_unifiedMedianPlain         = NULL;
   }
   if (m_subSampleVariancePlain) {
     delete m_subSampleVariancePlain;
     m_subSampleVariancePlain     = NULL;
   }
   if (m_unifiedSampleVariancePlain) {
     delete m_unifiedSampleVariancePlain;
     m_unifiedSampleVariancePlain = NULL;
   }
 
   return;
 }

template<class T >

const BaseEnvironment & QUESO::ScalarSequence< T >::env ( ) const

Access to QUESO environment.

Definition at line 101 of file ScalarSequence.C.

Referenced by QUESO::ComputeCovCorrBetweenScalarSequences(), and QUESO::ComputeUnifiedGaussian2dKde().

 {
   return m_env;
 }

template<class T >

void QUESO::ScalarSequence< T >::erasePositions	(	unsigned int	initialPos,
		unsigned int	numPos
	)

Erases numPos values of the sequence, starting at position initialPos.

This routine deletes all stored computed scalars.

Definition at line 206 of file ScalarSequence.C.

References queso_require_equal_to_msg, and queso_require_msg.

Referenced by QUESO::ArrayOfSequences< V, M >::erasePositions().

 {
   if (this->subSequenceSize() == 0) return;
 
   bool bRC = ((initialPos          <  this->subSequenceSize()) &&
               (0                   <  numPos                 ) &&
               ((initialPos+numPos) <= this->subSequenceSize()));
   queso_require_msg(bRC, "invalid input data");
 
   seqScalarPositionIteratorTypedef posIteratorBegin = m_seq.begin();
   if (initialPos < this->subSequenceSize()) std::advance(posIteratorBegin,initialPos);
   else                                      posIteratorBegin = m_seq.end();
 
   unsigned int posEnd = initialPos + numPos - 1;
   seqScalarPositionIteratorTypedef posIteratorEnd = m_seq.begin();
   if (posEnd < this->subSequenceSize()) std::advance(posIteratorEnd,posEnd);
   else                                  posIteratorEnd = m_seq.end();
 
   unsigned int oldSequenceSize = this->subSequenceSize();
   m_seq.erase(posIteratorBegin,posIteratorEnd);
   queso_require_equal_to_msg((oldSequenceSize - numPos), this->subSequenceSize(), "(oldSequenceSize - numPos) != this->subSequenceSize()");
 
   deleteStoredScalars();
 
   return;
 }

template<class T >

void QUESO::ScalarSequence< T >::extractRawData	(	unsigned int	initialPos,
		unsigned int	spacing,
		unsigned int	numPos,
		std::vector< double > &	rawData
	)		const

private

Extracts the raw data.

This method saves in rawData the data from the sequence of scalars (in private attribute m_seq) starting at position (initialPos), with a spacing of spacing until numPos positions have been extracted.

Definition at line 3158 of file ScalarSequence.C.

 {
   rawDataVec.resize(numPos);
   if (spacing == 1) {
     for (unsigned int j = 0; j < numPos; ++j) {
       rawDataVec[j] = m_seq[initialPos+j        ];
     }
   }
   else {
     for (unsigned int j = 0; j < numPos; ++j) {
       rawDataVec[j] = m_seq[initialPos+j*spacing];
     }
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::extractScalarSeq	(	unsigned int	initialPos,
		unsigned int	spacing,
		unsigned int	numPos,
		ScalarSequence< T > &	scalarSeq
	)		const

private

Extracts a sequence of scalars.

The sequence of scalars has size numPos, and it will be extracted starting at position (initialPos) of this sequence of scalars, given spacing spacing.

Definition at line 3117 of file ScalarSequence.C.

References QUESO::ScalarSequence< T >::resizeSequence().

 {
   scalarSeq.resizeSequence(numPos);
   if (spacing == 1) {
     for (unsigned int j = 0; j < numPos; ++j) {
       //if ((initialPos+j*spacing) > m_seq.size()) {
       //  std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
       //            << ": initialPos = " << initialPos
       //            << ", spacing = "    << spacing
       //            << ", numPos = "     << numPos
       //            << ", j = "          << j
       //            << ", position got too large"
       //            << std::endl;
       //}
       scalarSeq[j] = m_seq[initialPos+j        ];
     }
   }
   else {
     for (unsigned int j = 0; j < numPos; ++j) {
       //if ((initialPos+j*spacing) > m_seq.size()) {
       //  std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
       //            << ": initialPos = " << initialPos
       //            << ", spacing = "    << spacing
       //            << ", numPos = "     << numPos
       //            << ", j = "          << j
       //            << ", position got too large"
       //            << std::endl;
       //}
       scalarSeq[j] = m_seq[initialPos+j*spacing];
     }
   }
 
   return;
 }

template<class T >

void QUESO::ScalarSequence< T >::filter	(	unsigned int	initialPos,
		unsigned int	spacing
	)

Filters positions in the sequence of vectors.

Filtered positions will start at initialPos, and with spacing given by spacing.

Definition at line 2429 of file ScalarSequence.C.

Referenced by QUESO::MLSampling< P_V, P_M >::generateSequence(), QUESO::MetropolisHastingsSG< P_V, P_M >::generateSequence(), and QUESO::MLSampling< P_V, P_M >::generateSequence_Step11_inter0().

 {
   if (m_env.subDisplayFile()) {
     *m_env.subDisplayFile() << "Entering ScalarSequence<V,M>::filter()"
                             << ": initialPos = "      << initialPos
                             << ", spacing = "         << spacing
                             << ", subSequenceSize = " << this->subSequenceSize()
                             << std::endl;
   }
 
   unsigned int i = 0;
   unsigned int j = initialPos;
   unsigned int originalSubSequenceSize = this->subSequenceSize();
   while (j < originalSubSequenceSize) {
     if (i != j) {
       //*m_env.subDisplayFile() << i << "--" << j << " ";
       m_seq[i] = m_seq[j];
     }
     i++;
     j += spacing;
   }
 
   this->resizeSequence(i);
 
   if (m_env.subDisplayFile()) {
     *m_env.subDisplayFile() << "Leaving ScalarSequence<V,M>::filter()"
                             << ": initialPos = "      << initialPos
                             << ", spacing = "         << spacing
                             << ", subSequenceSize = " << this->subSequenceSize()
                             << std::endl;
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::getUnifiedContentsAtProc0Only	(	bool	useOnlyInter0Comm,
		std::vector< T > &	outputVec
	)		const

Gets the unified contents of processor of rank equals to 0.

Definition at line 235 of file ScalarSequence.C.

References queso_error_msg, queso_require_equal_to_msg, RawValue_MPI_DOUBLE, RawValue_MPI_IN_PLACE, and RawValue_MPI_INT.

Referenced by QUESO::MLSampling< P_V, P_M >::generateSequence_Step05_inter0(), and QUESO::MLSampling< P_V, P_M >::generateSequence_Step09_all().

 {
   // The logic (numSubEnvs == 1) does *not* apply here because 'outputVec' needs to be filled
   //if (m_env.numSubEnvironments() == 1) {
   //  // No need to do anything
   //  return;
   //}
 
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       int auxSubSize = (int) this->subSequenceSize();
       unsigned int auxUnifiedSize = this->unifiedSequenceSize(useOnlyInter0Comm);
       outputVec.resize(auxUnifiedSize,0.);
 
       //******************************************************************
       // Use MPI_Gatherv for the case different nodes have different amount of data // KAUST4
       //******************************************************************
       std::vector<int> recvcnts(m_env.inter0Comm().NumProc(),0); // '0' is NOT the correct value for recvcnts[0]
       m_env.inter0Comm().Gather((void *) &auxSubSize, 1, RawValue_MPI_INT, (void *) &recvcnts[0], (int) 1, RawValue_MPI_INT, 0,
                                 "ScalarSequence<T>::getUnifiedContentsAtProc0Only()",
                                 "failed MPI.Gather()");
       if (m_env.inter0Rank() == 0) {
         //recvcnts[0] = (int) this->subSequenceSize(); // FIX ME: really necessary????
         queso_require_equal_to_msg(recvcnts[0], (int) this->subSequenceSize(), "failed MPI.Gather() result at proc 0");
       }
 
       std::vector<int> displs(m_env.inter0Comm().NumProc(),0);
       for (unsigned int r = 1; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) { // Yes, from '1' on
         displs[r] = displs[r-1] + recvcnts[r-1];
       }
 
 #if 0 // for debug only
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
         for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) {
           *m_env.subDisplayFile() << "  auxSubSize = "            << auxSubSize
                                   << ", recvcnts[" << r << "] = " << recvcnts[r]
                                   << ", displs["   << r << "] = " << displs[r]
                                   << ", m_seq.size() = "          << m_seq.size()
                                   << ", outputVec.size() = "      << outputVec.size()
                                   << std::endl;
         }
         for (unsigned int i = 0; i < m_seq.size(); ++i) {
           *m_env.subDisplayFile() << "  (before gatherv) m_seq[" << i << "]= " << m_seq[i]
                                   << std::endl;
         }
       }
 #endif
       m_env.inter0Comm().Gatherv((void *) &m_seq[0], auxSubSize, RawValue_MPI_DOUBLE, (void *) &outputVec[0], (int *) &recvcnts[0], (int *) &displs[0], RawValue_MPI_DOUBLE, 0,
                                  "ScalarSequence<T>::getUnifiedContentsAtProc0Only()",
                                  "failed MPI.Gatherv()");
 
 #if 0 // for debug only
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
         for (unsigned int i = 0; i < m_seq.size(); ++i) {
           *m_env.subDisplayFile() << "  (after gatherv) m_seq[" << i << "]= " << m_seq[i]
                                   << std::endl;
         }
         for (unsigned int i = 0; i < outputVec.size(); ++i) {
           *m_env.subDisplayFile() << "  (after gatherv) outputVec[" << i << "]= " << outputVec[i]
                                   << std::endl;
         }
       }
 #endif
 
 #if 0 // for debug only
       if (m_env.inter0Rank() == 0) {
         for (unsigned int i = 0; i < auxSubSize; ++i) {
           outputVec[i] = m_seq[i];
         }
         m_env.inter0Comm().Gatherv(RawValue_MPI_IN_PLACE, auxSubSize, RawValue_MPI_DOUBLE, (void *) &outputVec[0], (int *) &recvcnts[0], (int *) &displs[0], RawValue_MPI_DOUBLE, 0,
                                    "ScalarSequence<T>::getUnifiedContentsAtProc0Only(1)",
                                    "failed MPI.Gatherv()");
       }
       else {
         m_env.inter0Comm().Gatherv((void *) &m_seq[0], auxSubSize, RawValue_MPI_DOUBLE, (void *) &outputVec[0], (int *) &recvcnts[0], (int *) &displs[0], RawValue_MPI_DOUBLE, 0,
                                    "ScalarSequence<T>::getUnifiedContentsAtProc0Only(2)",
                                    "failed MPI.Gatherv()");
       }
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
         for (unsigned int i = 0; i < m_seq.size(); ++i) {
           *m_env.subDisplayFile() << "  (after 2nd gatherv) m_seq[" << i << "]= " << m_seq[i]
                                   << std::endl;
         }
         for (unsigned int i = 0; i < outputVec.size(); ++i) {
           *m_env.subDisplayFile() << "  (after 2nd gatherv) outputVec[" << i << "]= " << outputVec[i]
                                   << std::endl;
         }
       }
 #endif
     }
     else {
       // Node not in the 'inter0' communicator
       // No need to do anything
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   return;
 }

template<class T >

const std::string & QUESO::ScalarSequence< T >::name ( ) const

Access to the name of the sequence of scalars.

Definition at line 108 of file ScalarSequence.C.

 {
   return m_name;
 }

template<class T>

ScalarSequence< T > & QUESO::ScalarSequence< T >::operator= ( const ScalarSequence< T > & rhs )

Assignment operator; it copies rhs to this.

Definition at line 61 of file ScalarSequence.C.

References copy.

 {
   this->copy(rhs);
   return *this;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::operator[] ( unsigned int posId ) const

Access position posId of the sequence of scalars (const).

Definition at line 69 of file ScalarSequence.C.

References queso_require_less_msg.

 {
   if (posId >= this->subSequenceSize()) {
     std::cerr << "In ScalarSequence<T>::operator[]() const"
               << ": posId = "                   << posId
               << ", this->subSequenceSize() = " << this->subSequenceSize()
               << std::endl;
   }
   queso_require_less_msg(posId, this->subSequenceSize(), "posId > subSequenceSize()");
 
   return m_seq[posId];
 }

template<class T >

T & QUESO::ScalarSequence< T >::operator[] ( unsigned int posId )

Access position posId of the sequence of scalars (non-const).

This routine deletes all stored computed scalars.

Definition at line 84 of file ScalarSequence.C.

References queso_require_less_msg.

 {
   if (posId >= this->subSequenceSize()) {
     std::cerr << "In ScalarSequence<T>::operator[]()"
               << ": posId = "                   << posId
               << ", this->subSequenceSize() = " << this->subSequenceSize()
               << std::endl;
   }
   queso_require_less_msg(posId, this->subSequenceSize(), "posId > subSequenceSize()");
 
   deleteStoredScalars();
 
   return m_seq[posId];
 }

template<class T>

void QUESO::ScalarSequence< T >::parallelMerge	(	std::vector< T > &	sortedBuffer,
		const std::vector< T > &	leafData,
		unsigned int	treeLevel
	)		const

private

Sorts/merges data in parallel using MPI.

Definition at line 3199 of file ScalarSequence.C.

References k, RawValue_MPI_DOUBLE, RawValue_MPI_UNSIGNED, SCALAR_SEQUENCE_DATA_MPI_MSG, SCALAR_SEQUENCE_INIT_MPI_MSG, and SCALAR_SEQUENCE_SIZE_MPI_MSG.

 {
   int parentNode = m_env.inter0Rank() & ~(1 << currentTreeLevel);
 
   if (m_env.inter0Rank() >= 0) { // KAUST
   if (currentTreeLevel == 0) {
     // Leaf node: sort own local data.
     unsigned int leafDataSize = leafData.size();
     sortedBuffer.resize(leafDataSize,0.);
     for (unsigned int i = 0; i < leafDataSize; ++i) {
       sortedBuffer[i] = leafData[i];
     }
     std::sort(sortedBuffer.begin(), sortedBuffer.end());
     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
       *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
                               << ": tree node "                                            << m_env.inter0Rank()
                               << ", leaf sortedBuffer[0] = "                               << sortedBuffer[0]
                               << ", leaf sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
                               << std::endl;
     }
   }
   else {
     int nextTreeLevel  = currentTreeLevel - 1;
     int rightChildNode = m_env.inter0Rank() | (1 << nextTreeLevel);
 
     if (rightChildNode >= m_env.inter0Comm().NumProc()) { // No right child. Move down one level.
       this->parallelMerge(sortedBuffer,
                           leafData,
                           nextTreeLevel);
     }
     else {
       unsigned int uintBuffer[1];
       uintBuffer[0] = nextTreeLevel;
       m_env.inter0Comm().Send((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, rightChildNode, SCALAR_SEQUENCE_INIT_MPI_MSG,
                               "ScalarSequence<T>::parallelMerge()",
                               "failed MPI.Send() for init");
 
       this->parallelMerge(sortedBuffer,
                           leafData,
                           nextTreeLevel);
 
       // Prepare variable 'leftSortedBuffer': just copy own current sorted data.
       unsigned int leftSize = sortedBuffer.size();
       std::vector<T> leftSortedBuffer(leftSize,0.);
       for (unsigned int i = 0; i < leftSize; ++i) {
         leftSortedBuffer[i] = sortedBuffer[i];
       }
 
       // Prepare variable 'rightSortedBuffer': receive data from right child node.
       RawType_MPI_Status status;
       m_env.inter0Comm().Recv((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, rightChildNode, SCALAR_SEQUENCE_SIZE_MPI_MSG, &status,
                               "ScalarSequence<T>::parallelMerge()",
                               "failed MPI.Recv() for size");
       //if (status) {}; // just to remove compiler warning
 
       unsigned int rightSize = uintBuffer[0];
       std::vector<T> rightSortedBuffer(rightSize,0.);
       m_env.inter0Comm().Recv((void *) &rightSortedBuffer[0], (int) rightSize, RawValue_MPI_DOUBLE, rightChildNode, SCALAR_SEQUENCE_DATA_MPI_MSG, &status,
                               "ScalarSequence<T>::parallelMerge()",
                               "failed MPI.Recv() for data");
 
       // Merge the two results into 'sortedBuffer'.
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
                                 << ": tree node "        << m_env.inter0Rank()
                                 << " is combining "      << leftSortedBuffer.size()
                                 << " left doubles with " << rightSortedBuffer.size()
                                 << " right doubles"
                                 << std::endl;
       }
 
       sortedBuffer.clear();
       sortedBuffer.resize(leftSortedBuffer.size()+rightSortedBuffer.size(),0.);
       unsigned int i = 0;
       unsigned int j = 0;
       unsigned int k = 0;
       while ((i < leftSize ) &&
              (j < rightSize)) {
         if (leftSortedBuffer[i] > rightSortedBuffer[j]) sortedBuffer[k++] = rightSortedBuffer[j++];
         else                                            sortedBuffer[k++] = leftSortedBuffer [i++];
       }
       while (i < leftSize ) sortedBuffer[k++] = leftSortedBuffer [i++];
       while (j < rightSize) sortedBuffer[k++] = rightSortedBuffer[j++];
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
                                 << ": tree node "                                              << m_env.inter0Rank()
                                 << ", merged sortedBuffer[0] = "                               << sortedBuffer[0]
                                 << ", merged sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
                                 << std::endl;
       }
     }
   }
 
   if (parentNode != m_env.inter0Rank()) {
     // Transmit data to parent node.
     unsigned int uintBuffer[1];
     uintBuffer[0] = sortedBuffer.size();
     m_env.inter0Comm().Send((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, parentNode, SCALAR_SEQUENCE_SIZE_MPI_MSG,
                             "ScalarSequence<T>::parallelMerge()",
                             "failed MPI.Send() for size");
 
     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
       *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
                               << ": tree node "                                          << m_env.inter0Rank()
                               << " is sending "                                          << sortedBuffer.size()
                               << " doubles to tree node "                                << parentNode
                               << ", with sortedBuffer[0] = "                             << sortedBuffer[0]
                               << " and sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
                               << std::endl;
     }
 
     m_env.inter0Comm().Send((void *) &sortedBuffer[0], (int) sortedBuffer.size(), RawValue_MPI_DOUBLE, parentNode, SCALAR_SEQUENCE_DATA_MPI_MSG,
                             "ScalarSequence<T>::parallelMerge()",
                             "failed MPI.Send() for data");
   }
   } // KAUST
 
   return;
 }

template<class T >

std::vector< T > & QUESO::ScalarSequence< T >::rawData ( )

private

The sequence of scalars. Access to private attribute m_seq.

Definition at line 3181 of file ScalarSequence.C.

Referenced by QUESO::ScalarSequence< T >::unifiedSort().

 {
   return m_seq;
 }

template<class T >

void QUESO::ScalarSequence< T >::resetValues	(	unsigned int	initialPos,
		unsigned int	numPos
	)

Sets numPos values of the sequence to zero, starting at position initialPos.

This routine deletes all stored computed scalars.

Definition at line 186 of file ScalarSequence.C.

References queso_require_msg.

 {
   if (this->subSequenceSize() == 0) return;
 
   bool bRC = ((initialPos          <  this->subSequenceSize()) &&
               (0                   <  numPos                 ) &&
               ((initialPos+numPos) <= this->subSequenceSize()));
   queso_require_msg(bRC, "invalid input data");
 
   for (unsigned int j = 0; j < numPos; ++j) {
     m_seq[initialPos+j] = 0.;
   }
 
   deleteStoredScalars();
 
   return;
 }

template<class T >

void QUESO::ScalarSequence< T >::resizeSequence ( unsigned int newSequenceSize )

Resizes the size of the sequence of scalars.

This routine deletes all stored computed scalars.

Definition at line 173 of file ScalarSequence.C.

Referenced by QUESO::ArrayOfSequences< V, M >::extractScalarSeq(), QUESO::SequenceOfVectors< V, M >::extractScalarSeq(), QUESO::ScalarSequence< T >::extractScalarSeq(), QUESO::MetropolisHastingsSG< P_V, P_M >::generateFullChain(), QUESO::MLSampling< P_V, P_M >::generateSequence_Level0_all(), QUESO::ScalarSequence< T >::subMedianExtra(), QUESO::ScalarSequence< T >::subPositionsOfMaximum(), QUESO::ScalarSequence< T >::subSort(), QUESO::ScalarSequence< T >::unifiedPositionsOfMaximum(), and QUESO::ScalarSequence< T >::unifiedSort().

 {
   if (newSequenceSize != this->subSequenceSize()) {
     m_seq.resize(newSequenceSize,0.);
     std::vector<T>(m_seq).swap(m_seq);
     deleteStoredScalars();
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::setGaussian	(	const T &	mean,
		const T &	stdDev
	)

Sets the values of the sequence as a Gaussian distribution of mean given by meanVec and standard deviation by stdDevVec.

This routine deletes all stored computed scalars.

Definition at line 513 of file ScalarSequence.C.

Referenced by QUESO::ArrayOfSequences< V, M >::setGaussian().

 {
   unsigned int maxJ = this->subSequenceSize();
   if (meanValue == 0.) {
     for (unsigned int j = 0; j < maxJ; ++j) {
       m_seq[j] = m_env.rngObject()->gaussianSample(stdDev);
     }
   }
   else {
     for (unsigned int j = 0; j < maxJ; ++j) {
       m_seq[j] = meanValue + m_env.rngObject()->gaussianSample(stdDev);
     }
   }
 
   deleteStoredScalars();
 
   return;
 }

template<class T >

void QUESO::ScalarSequence< T >::setName ( const std::string & newName )

Sets a new name to the sequence of scalars.

Definition at line 115 of file ScalarSequence.C.

Referenced by QUESO::MLSampling< P_V, P_M >::generateSequence(), QUESO::MLSampling< P_V, P_M >::generateSequence_Level0_all(), QUESO::MLSampling< P_V, P_M >::generateSequence_Step02_inter0(), and QUESO::MLSampling< P_V, P_M >::generateSequence_Step11_inter0().

 {
   m_name = newName;
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::setUniform	(	const T &	a,
		const T &	b
	)

Sets the values of the sequence as a uniform distribution between the values given by vectors aVec and bVec.

This routine deletes all stored computed scalars.

Definition at line 534 of file ScalarSequence.C.

Referenced by QUESO::ArrayOfSequences< V, M >::setUniform().

 {
   unsigned int maxJ = this->subSequenceSize();
   if (a == 0.) {
     if (b == 1.) {
       for (unsigned int j = 0; j < maxJ; ++j) {
         m_seq[j] = m_env.rngObject()->uniformSample();
       }
     }
     else {
       for (unsigned int j = 0; j < maxJ; ++j) {
         m_seq[j] = b*m_env.rngObject()->uniformSample();
       }
     }
   }
   else {
     if ((b-a) == 1.) {
       for (unsigned int j = 0; j < maxJ; ++j) {
         m_seq[j] = a + m_env.rngObject()->uniformSample();
       }
     }
     else {
       for (unsigned int j = 0; j < maxJ; ++j) {
         m_seq[j] = a + (b-a)*m_env.rngObject()->uniformSample();
       }
     }
   }
 
   deleteStoredScalars();
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subBasicCdf	(	unsigned int	numIntervals,
		UniformOneDGrid< T > *&	gridValues,
		std::vector< T > &	cdfValues
	)		const

Finds the Cumulative Distribution Function (CDF) of the sub-sequence of scalars.

Definition at line 710 of file ScalarSequence.C.

 {
   T                         tmpMinValue;
   T                         tmpMaxValue;
   std::vector<unsigned int> bins(numEvaluationPoints,0);
 
   subMinMaxExtra(0, // initialPos
                  this->subSequenceSize(),
                  tmpMinValue,
                  tmpMaxValue);
   subBasicHistogram(0, // initialPos,
                     tmpMinValue,
                     tmpMaxValue,
                     gridValues,
                     bins);
 
   unsigned int sumOfBins = 0;
   for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
     sumOfBins += bins[i];
   }
 
   cdfValues.clear();
   cdfValues.resize(numEvaluationPoints);
   unsigned int partialSum = 0;
   for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
     partialSum += bins[i];
     cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subBasicHistogram	(	unsigned int	initialPos,
		const T &	minHorizontalValue,
		const T &	maxHorizontalValue,
		UniformOneDGrid< T > *&	gridValues,
		std::vector< unsigned int > &	bins
	)		const

Calculates the histogram of the sub-sequence.

It requires the specification of the maximum (maxHorizontalValue) and the minimum (minHorizontalValue) values if the data and the initial position (initialPos) from where the data will be considered. Output: grid values that will act as the center of each bin (gridValues) and the number of data occurrences in each bin (bins).

Definition at line 1747 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

 {
   queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
 
   for (unsigned int j = 0; j < bins.size(); ++j) {
     bins[j] = 0;
   }
 
   double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
   double minCenter = minHorizontalValue - horizontalDelta/2.;
   double maxCenter = maxHorizontalValue + horizontalDelta/2.;
   gridValues = new UniformOneDGrid<T>(m_env,
                                              "",
                                              bins.size(),
                                        minCenter,
                                        maxCenter);
 
   unsigned int dataSize = this->subSequenceSize();
   for (unsigned int j = 0; j < dataSize; ++j) {
     double value = m_seq[j];
     if (value < minHorizontalValue) {
       bins[0] += value;
     }
     else if (value >= maxHorizontalValue) {
       bins[bins.size()-1] += value;
     }
     else {
       unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
       bins[index] += value;
     }
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subGaussian1dKde	(	unsigned int	initialPos,
		double	scaleValue,
		const std::vector< T > &	evaluationPositions,
		std::vector< double > &	densityValues
	)		const

Gaussian kernel for the KDE estimate of the sub-sequence.

Computes a probability density estimate of the sample in this sub-sequence, starting at position initialPos. densityValues is the vector of density values evaluated at the points in evaluationPositions. The estimate is based on Gaussian (normal) kernel function, using a window parameter (scaleValue).

Definition at line 2319 of file ScalarSequence.C.

References k, QUESO::MiscGaussianDensity(), and queso_require_msg.

Referenced by QUESO::SequenceOfVectors< V, M >::subGaussian1dKde().

 {
   bool bRC = ((initialPos                 <  this->subSequenceSize()   ) &&
               (0                          <  evaluationPositions.size()) &&
               (evaluationPositions.size() == densityValues.size()      ));
   queso_require_msg(bRC, "invalid input data");
 
   unsigned int dataSize = this->subSequenceSize() - initialPos;
   unsigned int numEvals = evaluationPositions.size();
 
   double scaleInv = 1./scaleValue;
   for (unsigned int j = 0; j < numEvals; ++j) {
     double x = evaluationPositions[j];
     double value = 0.;
     for (unsigned int k = 0; k < dataSize; ++k) {
       double xk = m_seq[initialPos+k];
       value += MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
     }
     densityValues[j] = scaleInv * (value/(double) dataSize);
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subHistogram	(	unsigned int	initialPos,
		const T &	minHorizontalValue,
		const T &	maxHorizontalValue,
		std::vector< T > &	centers,
		std::vector< unsigned int > &	bins
	)		const

Calculates the histogram of the sub-sequence.

It requires the specification of the maximum (maxHorizontalValue) and the minimum (minHorizontalValue) values if the data and the initial position (initialPos) from where the data will be considered. Output: the center of each bin (centers) and the number of data occurrences in each bin (bins).

Definition at line 1610 of file ScalarSequence.C.

References queso_require_equal_to_msg, and queso_require_greater_equal_msg.

Referenced by QUESO::SequenceOfVectors< V, M >::subHistogram().

 {
   queso_require_equal_to_msg(centers.size(), bins.size(), "vectors 'centers' and 'bins' have different sizes");
 
   queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
 
   if (initialPos) {}; // just to remove compiler warning
 
   for (unsigned int j = 0; j < bins.size(); ++j) {
     centers[j] = 0.;
     bins[j] = 0;
   }
 
   double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
 
   double minCenter = minHorizontalValue - horizontalDelta/2.;
   double maxCenter = maxHorizontalValue + horizontalDelta/2.;
   for (unsigned int j = 0; j < centers.size(); ++j) {
     double factor = ((double) j)/(((double) centers.size()) - 1.);
     centers[j] = (1. - factor) * minCenter + factor * maxCenter;
   }
 
   unsigned int dataSize = this->subSequenceSize();
   for (unsigned int j = 0; j < dataSize; ++j) {
     double value = m_seq[j];
     if (value < minHorizontalValue) {
       bins[0]++;
     }
     else if (value >= maxHorizontalValue) {
       bins[bins.size()-1]++;
     }
     else {
       unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
       bins[index]++;
     }
   }
 
   return;
 }

template<class T >

T QUESO::ScalarSequence< T >::subInterQuantileRange ( unsigned int initialPos ) const

Returns the interquartile range of the values in the sub-sequence.

The IQR is a robust estimate of the spread of the data, since changes in the upper and lower 25% of the data do not affect it. If there are outliers in the data, then the IQR is more representative than the standard deviation as an estimate of the spread of the body of the data. The IQR is less efficient than the standard deviation as an estimate of the spread when the data is all from the normal distribution. (from Matlab)

Definition at line 2004 of file ScalarSequence.C.

References queso_require_equal_to_msg, queso_require_less_msg, and QUESO::ScalarSequence< T >::subSequenceSize().

Referenced by QUESO::SequenceOfVectors< V, M >::subInterQuantileRange().

 {
   queso_require_less_msg(initialPos, this->subSequenceSize(), "'initialPos' is too big");
 
   ScalarSequence sortedSequence(m_env,0,"");
   this->subSort(initialPos,
                 sortedSequence);
 
   // The test above guarantees that 'dataSize >= 1'
   unsigned int dataSize = this->subSequenceSize() - initialPos;
 
   queso_require_equal_to_msg(dataSize, sortedSequence.subSequenceSize(), "inconsistent size variables");
 
   bool everythingOk = true;
 
   // pos1 = (dataSize+1)/4 - 1
   // pos1 >= 0            <==> dataSize   >= 3
   // pos1 <  (dataSize-1) <==> 3*dataSize >  1
   unsigned int pos1 = (unsigned int) ( (((double) dataSize) + 1.)*1./4. - 1. );
   if (pos1 > (dataSize-1)) {
     pos1 = 0;
     everythingOk = false;
   }
   unsigned int pos1inc = pos1+1;
   if (pos1inc > (dataSize-1)) {
     pos1inc = dataSize-1;
     everythingOk = false;
   }
 
   // pos3 = (dataSize+1)*3/4 - 1
   // pos3 >= 0            <==> dataSize >= 1/3
   // pos3 <  (dataSize-1) <==> dataSize >  3
   unsigned int pos3 = (unsigned int) ( (((double) dataSize) + 1.)*3./4. - 1. );
   if (pos3 > (dataSize-1)) {
     pos3 = 0;
     everythingOk = false;
   }
   unsigned int pos3inc = pos3+1;
   if (pos3inc > (dataSize-1)) {
     pos3inc = dataSize-1;
     everythingOk = false;
   }
 
   double fraction1 = (((double) dataSize) + 1.)*1./4. - 1. - ((double) pos1);
   if (fraction1 < 0.) {
     fraction1 = 0.;
     everythingOk = false;
   }
   double fraction3 = (((double) dataSize) + 1.)*3./4. - 1. - ((double) pos3);
   if (fraction3 < 0.) {
     fraction3 = 0.;
     everythingOk = false;
   }
 
   if (everythingOk == false) {
     std::cerr << "In ScalarSequence<T>::subInterQuantileRange()"
               << ", worldRank = " << m_env.worldRank()
               << ": at least one adjustment was necessary"
               << std::endl;
   }
 
   //if (m_env.subDisplayFile()) {
   //  *m_env.subDisplayFile() << "In ScalarSequence::subInterQuantileRange()"
   //                          << ", initialPos = "               << initialPos
   //                          << ", this->subSequenceSize() = "  << this->subSequenceSize()
   //                          << ", dataSize = "                 << dataSize
   //                          << ", sortedSequence.size() = "    << sortedSequence.size()
   //                          << ", pos1 = "                     << pos1
   //                          << ", pos3 = "                     << pos3
   //                          << std::endl;
   //}
 
   T value1 = (1.-fraction1) * sortedSequence[pos1] + fraction1 * sortedSequence[pos1inc];
   T value3 = (1.-fraction3) * sortedSequence[pos3] + fraction3 * sortedSequence[pos3inc];
   T iqrValue = value3 - value1;
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
     *m_env.subDisplayFile() << "In ScalarSequence<T>::subInterQuantileRange()"
                             << ": iqrValue = " << iqrValue
                             << ", dataSize = " << dataSize
                             << ", pos1 = "     << pos1
                             << ", pos3 = "     << pos3
                             << ", value1 = "   << value1
                             << ", value3 = "   << value3
                             << std::endl;
 
     // Save data only once into a separate file
     //std::ofstream* ofsvar = new std::ofstream(("sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
     //if ((ofsvar            == NULL ) ||
     //    (ofsvar->is_open() == false)) {
     //  delete ofsvar;
     //  ofsvar = new std::ofstream(("sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
 
     //  *ofsvar << "var_sort_sub" << m_env.subIdString() << " = zeros(" << 1
     //          << ","                                                  << dataSize
     //          << ");"
     //          << std::endl;
     //  for (unsigned int j = 0; j < dataSize; ++j) {
     //    *ofsvar << "var_sort_sub" << m_env.subIdString() << "(" << 1
     //            << ","                                          << j+1
     //            << ") = "                                       << sortedSequence[j]
     //            << ";"
     //            << std::endl;
     //  }
     //}
     //delete ofsvar;
   }
 
   return iqrValue;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::subMaxPlain ( ) const

Finds the maximum value of the sub-sequence of scalars.

Definition at line 367 of file ScalarSequence.C.

Referenced by QUESO::BaseVectorSequence< V, M >::subPositionsOfMaximum(), QUESO::ScalarSequence< T >::subPositionsOfMaximum(), QUESO::BaseVectorSequence< V, M >::unifiedPositionsOfMaximum(), and QUESO::ScalarSequence< T >::unifiedPositionsOfMaximum().

 {
   if (m_subMaxPlain == NULL) {
     if (m_subMinPlain == NULL) m_subMinPlain = new T(0.);
     m_subMaxPlain = new T(0.);
     subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
   }
 
   return *m_subMaxPlain;
 }

template<class T >

T QUESO::ScalarSequence< T >::subMeanExtra	(	unsigned int	initialPos,
		unsigned int	numPos
	)		const

Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPos.

Definition at line 836 of file ScalarSequence.C.

References queso_require_msg.

Referenced by QUESO::ArrayOfSequences< V, M >::mean(), and QUESO::SequenceOfVectors< V, M >::subMeanExtra().

 {
   if (this->subSequenceSize() == 0) return 0.;
 
   bool bRC = ((initialPos          <  this->subSequenceSize()) &&
               (0                   <  numPos                 ) &&
               ((initialPos+numPos) <= this->subSequenceSize()));
   if (bRC == false) {
     std::cerr << "In ScalarSequence<T>::subMeanExtra()"
               << ": ERROR at fullRank "         << m_env.fullRank()
               << ", initialPos = "              << initialPos
               << ", numPos = "                  << numPos
               << ", this->subSequenceSize() = " << this->subSequenceSize()
               << std::endl;
     if (m_env.subDisplayFile()) {
       *m_env.subDisplayFile() << "In ScalarSequence<T>::subMeanExtra()"
                               << ": ERROR at fullRank "         << m_env.fullRank()
                               << ", initialPos = "              << initialPos
                               << ", numPos = "                  << numPos
                               << ", this->subSequenceSize() = " << this->subSequenceSize()
                               << std::endl;
     }
   }
   queso_require_msg(bRC, "invalid input data");
 
   unsigned int finalPosPlus1 = initialPos + numPos;
   T tmpSum = 0.;
   for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
     tmpSum += m_seq[j];
   }
 
   return tmpSum/(T) numPos;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::subMeanPlain ( ) const

Finds the mean value of the sub-sequence of scalars.

Definition at line 393 of file ScalarSequence.C.

 {
   if (m_subMeanPlain == NULL) {
     m_subMeanPlain = new T(0.);
     *m_subMeanPlain = subMeanExtra(0,subSequenceSize());
   }
 
   return *m_subMeanPlain;
 }

template<class T >

T QUESO::ScalarSequence< T >::subMedianExtra	(	unsigned int	initialPos,
		unsigned int	numPos
	)		const

Finds the median value of the sub-sequence, considering numPos positions starting at position initialPos.

Definition at line 952 of file ScalarSequence.C.

References queso_require_msg, QUESO::ScalarSequence< T >::resizeSequence(), and QUESO::ScalarSequence< T >::subSort().

Referenced by QUESO::SequenceOfVectors< V, M >::subMedianExtra().

 {
   if (this->subSequenceSize() == 0) return 0.;
 
   bool bRC = ((initialPos          <  this->subSequenceSize()) &&
               (0                   <  numPos                 ) &&
               ((initialPos+numPos) <= this->subSequenceSize()));
   if (bRC == false) {
     std::cerr << "In ScalarSequence<T>::subMedianExtra()"
               << ": ERROR at fullRank "         << m_env.fullRank()
               << ", initialPos = "              << initialPos
               << ", numPos = "                  << numPos
               << ", this->subSequenceSize() = " << this->subSequenceSize()
               << std::endl;
     if (m_env.subDisplayFile()) {
       *m_env.subDisplayFile() << "In ScalarSequence<T>::subMedianExtra()"
                               << ": ERROR at fullRank "         << m_env.fullRank()
                               << ", initialPos = "              << initialPos
                               << ", numPos = "                  << numPos
                               << ", this->subSequenceSize() = " << this->subSequenceSize()
                               << std::endl;
     }
   }
   queso_require_msg(bRC, "invalid input data");
 
   ScalarSequence sortedSequence(m_env,0,"");
   sortedSequence.resizeSequence(numPos);
   this->extractScalarSeq(initialPos,
                          1,
                          numPos,
                          sortedSequence);
   sortedSequence.subSort();
 
   unsigned int tmpPos = (unsigned int) (0.5 * (double) numPos);
   T resultValue = sortedSequence[tmpPos];
 
   return resultValue;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::subMedianPlain ( ) const

Finds the median value of the sub-sequence of scalars.

Definition at line 417 of file ScalarSequence.C.

 {
   if (m_subMedianPlain == NULL) {
     m_subMedianPlain = new T(0.);
     *m_subMedianPlain = subMedianExtra(0,subSequenceSize());
   }
 
   return *m_subMedianPlain;
 }

template<class T>

void QUESO::ScalarSequence< T >::subMinMaxExtra	(	unsigned int	initialPos,
		unsigned int	numPos,
		T &	minValue,
		T &	maxValue
	)		const

Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting at position initialPos.

Definition at line 1510 of file ScalarSequence.C.

References queso_require_less_equal_msg, and queso_require_msg.

Referenced by QUESO::MLSampling< P_V, P_M >::generateSequence(), QUESO::ArrayOfSequences< V, M >::minMax(), and QUESO::SequenceOfVectors< V, M >::subMinMaxExtra().

 {
   queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid input");
 
   seqScalarPositionConstIteratorTypedef pos1 = m_seq.begin();
   std::advance(pos1,initialPos);
 
   seqScalarPositionConstIteratorTypedef pos2 = m_seq.begin();
   std::advance(pos2,initialPos+numPos);
 
   if ((initialPos+numPos) == this->subSequenceSize()) {
     queso_require_msg(!(pos2 != m_seq.end()), "invalid state");
   }
 
   seqScalarPositionConstIteratorTypedef pos;
   pos = std::min_element(pos1, pos2);
   minValue = *pos;
   pos = std::max_element(pos1, pos2);
   maxValue = *pos;
 
   return;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::subMinPlain ( ) const

Finds the minimum value of the sub-sequence of scalars.

Definition at line 341 of file ScalarSequence.C.

 {
   if (m_subMinPlain == NULL) {
     m_subMinPlain = new T(0.);
     if (m_subMaxPlain == NULL) m_subMaxPlain = new T(0.);
     subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
   }
 
   return *m_subMinPlain;
 }

template<class T>

T QUESO::ScalarSequence< T >::subPopulationVariance	(	unsigned int	initialPos,
		unsigned int	numPos,
		const T &	meanValue
	)		const

Finds the population variance of the sub-sequence, considering numPos positions starting at position initialPos and of mean meanValue.

Output: calculated population variance of the sub-sequence of scalars. The population variance $\sigma^2$ is defined by $\sigma^2 = \frac{1}{n-1} \sum_{i=1}^n \left(y_i - \mu \right)^2$ , where $\mu$ is the sample mean and $ n $ is the sample size. This procedure lets the users choose the initial position and the number of elements of the sequence which will be used to evaluate the population variance .

Definition at line 1216 of file ScalarSequence.C.

References queso_require_msg.

Referenced by QUESO::SequenceOfVectors< V, M >::subPopulationVariance().

 {
   if (this->subSequenceSize() == 0) return 0.;
 
   bool bRC = ((initialPos          <  this->subSequenceSize()) &&
               (0                   <  numPos                 ) &&
               ((initialPos+numPos) <= this->subSequenceSize()));
   queso_require_msg(bRC, "invalid input data");
 
   unsigned int finalPosPlus1 = initialPos + numPos;
   T diff;
   T popValue = 0.;
   for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
     diff = m_seq[j] - meanValue;
     popValue += diff*diff;
   }
 
   popValue /= (T) numPos;
 
   return popValue;
 }

template<class T>

T QUESO::ScalarSequence< T >::subPositionsOfMaximum	(	const ScalarSequence< T > &	subCorrespondingScalarValues,
		ScalarSequence< T > &	subPositionsOfMaximum
	)

Finds the positions where the maximum element occurs in the sub-sequence.

Definition at line 2518 of file ScalarSequence.C.

References queso_require_equal_to_msg, QUESO::ScalarSequence< T >::resizeSequence(), QUESO::ScalarSequence< T >::subMaxPlain(), and QUESO::ScalarSequence< T >::subSequenceSize().

 {
   queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
 
   T subMaxValue = subCorrespondingScalarValues.subMaxPlain();
   unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
 
   unsigned int subNumPos = 0;
   for (unsigned int i = 0; i < iMax; ++i) {
     if (subCorrespondingScalarValues[i] == subMaxValue) {
       subNumPos++;
     }
   }
 
   subPositionsOfMaximum.resizeSequence(subNumPos);
   unsigned int j = 0;
   for (unsigned int i = 0; i < iMax; ++i) {
     if (subCorrespondingScalarValues[i] == subMaxValue) {
       subPositionsOfMaximum[j] = (*this)[i];
       j++;
     }
   }
 
   return subMaxValue;
 }

template<class T>

T QUESO::ScalarSequence< T >::subSampleStd	(	unsigned int	initialPos,
		unsigned int	numPos,
		const T &	meanValue
	)		const

Finds the sample standard deviation of the unified sequence, considering numPos positions starting at position initialPos and of mean meanValue.

Definition at line 1128 of file ScalarSequence.C.

References queso_require_msg.

Referenced by QUESO::SequenceOfVectors< V, M >::subSampleStd().

 {
   if (this->subSequenceSize() == 0) return 0.;
 
   bool bRC = ((initialPos          <  this->subSequenceSize()) &&
               (0                   <  numPos                 ) &&
               ((initialPos+numPos) <= this->subSequenceSize()));
   queso_require_msg(bRC, "invalid input data");
 
   unsigned int finalPosPlus1 = initialPos + numPos;
   T diff;
   T stdValue = 0.;
   for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
     diff = m_seq[j] - meanValue;
     stdValue += diff*diff;
   }
 
   stdValue /= (((T) numPos) - 1.);
   stdValue = sqrt(stdValue);
 
   return stdValue;
 }

template<class T>

T QUESO::ScalarSequence< T >::subSampleVarianceExtra	(	unsigned int	initialPos,
		unsigned int	numPos,
		const T &	meanValue
	)		const

Finds the sample variance of the sub-sequence, considering numPos positions starting at position initialPos and of mean meanVec.

Output: calculated sample variance of the sub-sequence of scalars. The sample variance $\sigma_y^2$ is the second sample central moment and is defined by $\sigma_y^2 = \frac{1}{n} \sum_{i=1}^n \left(y_i - \mu \right)^2$ , where $\mu$ is the sample mean and $ n $ is the sample size. This procedure lets the users choose the initial position and the number of elements of the sequence which will be used to evaluate the sample variance.

Definition at line 1042 of file ScalarSequence.C.

References queso_require_msg.

Referenced by QUESO::SequenceOfVectors< V, M >::subSampleVarianceExtra().

 {
   if (this->subSequenceSize() == 0) return 0.;
 
   bool bRC = ((initialPos          <  this->subSequenceSize()) &&
               (0                   <  numPos                 ) &&
               ((initialPos+numPos) <= this->subSequenceSize()));
   queso_require_msg(bRC, "invalid input data");
 
   unsigned int finalPosPlus1 = initialPos + numPos;
   T diff;
   T samValue = 0.;
   for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
     diff = m_seq[j] - meanValue;
     samValue += diff*diff;
   }
 
   samValue /= (((T) numPos) - 1.);
 
   return samValue;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::subSampleVariancePlain ( ) const

Finds the variance of a sample of the sub-sequence of scalars.

Definition at line 441 of file ScalarSequence.C.

 {
   if (m_subSampleVariancePlain == NULL) {
     m_subSampleVariancePlain = new T(0.);
     *m_subSampleVariancePlain = subSampleVarianceExtra(0,subSequenceSize(),subMeanPlain());
   }
 
   return *m_subSampleVariancePlain;
 }

template<class T>

T QUESO::ScalarSequence< T >::subScaleForKde	(	unsigned int	initialPos,
		const T &	iqrValue,
		unsigned int	kdeDimension
	)		const

Selects the scales (output value) for the kernel density estimation, considering only the sub-sequence.

The bandwidth of the kernel is a free parameter which exhibits a strong influence on the resulting estimate. Silverman (1986) suggests the following normal-based estimates: S1 = 1.06 × (standard deviation) × n^{-1/5} S2 = 0.79 × (iqrValue) × n^{-1/5}, where iqrValue is the interquartile range OutputValue = 0.90 × minimum(standard deviation, iqrValue /1.34) × n^{-1/5}. These estimates are popular due to their simplicity, and are used in QUESO with the adaptation of the exponent oven the sample size n (-1/5) with -1/(4 + kdeDimension) where kdeDimension is the KDE dimension.

Definition at line 2208 of file ScalarSequence.C.

References queso_require_msg.

Referenced by QUESO::SequenceOfVectors< V, M >::subScalesForKde().

 {
   bool bRC = (initialPos < this->subSequenceSize());
   queso_require_msg(bRC, "invalid input data");
 
   unsigned int dataSize = this->subSequenceSize() - initialPos;
 
   T meanValue = this->subMeanExtra(initialPos,
                                    dataSize);
 
   T samValue = this->subSampleVarianceExtra(initialPos,
                                             dataSize,
                                             meanValue);
 
   T scaleValue;
   if (iqrValue <= 0.) {
     scaleValue = 1.06*std::sqrt(samValue)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
    }
   else {
     scaleValue = 1.06*std::min(std::sqrt(samValue),iqrValue/1.34)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
   }
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
     *m_env.subDisplayFile() << "In ScalarSequence<T>::subScaleForKde()"
                             << ": iqrValue = "   << iqrValue
                             << ", meanValue = "  << meanValue
                             << ", samValue = "   << samValue
                             << ", dataSize = "   << dataSize
                             << ", scaleValue = " << scaleValue
                             << std::endl;
   }
 
   return scaleValue;
 }

template<class T >

unsigned int QUESO::ScalarSequence< T >::subSequenceSize ( ) const

Size of the sub-sequence of scalars.

Definition at line 136 of file ScalarSequence.C.

Referenced by QUESO::ScalarSequence< T >::append(), QUESO::ComputeCovCorrBetweenScalarSequences(), QUESO::ComputeSubGaussian2dKde(), QUESO::ScalarSequence< T >::copy(), QUESO::MLSampling< P_V, P_M >::generateBalLinkedChains_all(), QUESO::MLSampling< P_V, P_M >::generateSequence(), QUESO::MLSampling< P_V, P_M >::generateSequence_Step02_inter0(), QUESO::MLSampling< P_V, P_M >::generateSequence_Step03_inter0(), QUESO::MLSampling< P_V, P_M >::generateSequence_Step04_inter0(), QUESO::MLSampling< P_V, P_M >::generateSequence_Step05_inter0(), QUESO::MLSampling< P_V, P_M >::generateSequence_Step09_all(), QUESO::MLSampling< P_V, P_M >::generateUnbLinkedChains_all(), QUESO::SubF1F2Gaussian2dKdeIntegral(), QUESO::ScalarSequence< T >::subInterQuantileRange(), QUESO::BaseVectorSequence< V, M >::subPositionsOfMaximum(), QUESO::ScalarSequence< T >::subPositionsOfMaximum(), QUESO::ScalarSequence< T >::unifiedInterQuantileRange(), QUESO::BaseVectorSequence< V, M >::unifiedPositionsOfMaximum(), QUESO::ScalarSequence< T >::unifiedPositionsOfMaximum(), and QUESO::ScalarSequence< T >::unifiedSort().

 {
   return m_seq.size();
 }

template<class T>

void QUESO::ScalarSequence< T >::subSort	(	unsigned int	initialPos,
		ScalarSequence< T > &	sortedSequence
	)		const

Sorts the sub-sequence of scalars.

Definition at line 1883 of file ScalarSequence.C.

References QUESO::ScalarSequence< T >::resizeSequence(), and QUESO::ScalarSequence< T >::subSort().

Referenced by QUESO::ScalarSequence< T >::subMedianExtra(), and QUESO::ScalarSequence< T >::subSort().

 {
   unsigned int numPos = this->subSequenceSize() - initialPos;
   sortedSequence.resizeSequence(numPos);
   this->extractScalarSeq(initialPos,
                          1,
                          numPos,
                          sortedSequence);
   sortedSequence.subSort();
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subSort ( )

private

Sorts the sequence of scalars in the private attribute m_seq.

Definition at line 3189 of file ScalarSequence.C.

 {
   std::sort(m_seq.begin(), m_seq.end());
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subUniformlySampledCdf	(	unsigned int	numIntervals,
		T &	minDomainValue,
		T &	maxDomainValue,
		std::vector< T > &	cdfValues
	)		const

Uniformly samples from the CDF from the sub-sequence.

Definition at line 569 of file ScalarSequence.C.

 {
   T                         tmpMinValue;
   T                         tmpMaxValue;
   std::vector<T>            centers(numEvaluationPoints,0.);
   std::vector<unsigned int> bins   (numEvaluationPoints,0);
 
   subMinMaxExtra(0, // initialPos
                  this->subSequenceSize(),
                  tmpMinValue,
                  tmpMaxValue);
   subHistogram(0, // initialPos,
                tmpMinValue,
                tmpMaxValue,
                centers,
                bins);
 
   minDomainValue = centers[0];
   maxDomainValue = centers[centers.size()-1];
 
   unsigned int sumOfBins = 0;
   for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
     sumOfBins += bins[i];
   }
 
   cdfValues.clear();
   cdfValues.resize(numEvaluationPoints);
   unsigned int partialSum = 0;
   for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
     partialSum += bins[i];
     cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subWeightCdf	(	unsigned int	numIntervals,
		std::vector< T > &	gridValues,
		std::vector< T > &	cdfValues
	)		const

Finds the Weighted Cumulative Distribution Function (CDF) of the sub-sequence of scalars.

Definition at line 747 of file ScalarSequence.C.

 {
   T                         tmpMinValue;
   T                         tmpMaxValue;
   std::vector<unsigned int> bins(numEvaluationPoints,0);
   gridValues.resize             (numEvaluationPoints,0.);
   cdfValues.resize              (numEvaluationPoints,0.);
 
   subMinMaxExtra(0, // initialPos
                  this->subSequenceSize(),
                  tmpMinValue,
                  tmpMaxValue);
 
   if (tmpMinValue == tmpMaxValue) {
     if (tmpMinValue < -1.e-12) {
       tmpMinValue += tmpMinValue*(1.e-8);
     }
     else if (tmpMinValue > 1.e-12) {
       tmpMinValue -= tmpMinValue*(1.e-8);
     }
     else {
       tmpMinValue = 1.e-12;
     }
   }
 
   subWeightHistogram(0, // initialPos,
                      tmpMinValue,
                      tmpMaxValue,
                      gridValues,
                      bins);
 
   unsigned int sumOfBins = 0;
   for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
     sumOfBins += bins[i];
   }
 
   cdfValues.clear();
   cdfValues.resize(numEvaluationPoints);
   unsigned int partialSum = 0;
   for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
     partialSum += bins[i];
     cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subWeightCdf	(	unsigned int	numIntervals,
		UniformOneDGrid< T > *&	gridValues,
		std::vector< T > &	cdfValues
	)		const

Finds the Weighted Cumulative Distribution Function (CDF) of the sub-sequence of scalars.

Definition at line 799 of file ScalarSequence.C.

 {
   T                         tmpMinValue;
   T                         tmpMaxValue;
   std::vector<unsigned int> bins(numEvaluationPoints,0);
 
   subMinMaxExtra(0, // initialPos
                  this->subSequenceSize(),
                  tmpMinValue,
                  tmpMaxValue);
   subWeightHistogram(0, // initialPos,
                      tmpMinValue,
                      tmpMaxValue,
                      gridValues,
                      bins);
 
   unsigned int sumOfBins = 0;
   for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
     sumOfBins += bins[i];
   }
 
   cdfValues.clear();
   cdfValues.resize(numEvaluationPoints);
   unsigned int partialSum = 0;
   for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
     partialSum += bins[i];
     cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subWeightHistogram	(	unsigned int	initialPos,
		const T &	minHorizontalValue,
		const T &	maxHorizontalValue,
		UniformOneDGrid< T > *&	gridValues,
		std::vector< unsigned int > &	bins
	)		const

Calculates the weighted histogram of the sub-sequence.

TODO:Note that this method is exactly the same as subBasicHistogram().

It requires the specification of the maximum (maxHorizontalValue) and the minimum (minHorizontalValue) values if the data and the initial position (initialPos) from where the data will be considered. Output: grid values that will act as the center of each bin (gridValues) and the number of data occurrences in each bin (bins).

Definition at line 1789 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

 {
   queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
 
   for (unsigned int j = 0; j < bins.size(); ++j) {
     bins[j] = 0;
   }
 
   double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
   double minCenter = minHorizontalValue - horizontalDelta/2.;
   double maxCenter = maxHorizontalValue + horizontalDelta/2.;
   gridValues = new UniformOneDGrid<T>(m_env,
                                              "",
                                              bins.size(),
                                        minCenter,
                                        maxCenter);
 
   unsigned int dataSize = this->subSequenceSize();
   for (unsigned int j = 0; j < dataSize; ++j) {
     double value = m_seq[j];
     if (value < minHorizontalValue) {
       bins[0] += value;
     }
     else if (value >= maxHorizontalValue) {
       bins[bins.size()-1] += value;
     }
     else {
       unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
       bins[index] += value;
     }
   }
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::subWeightHistogram	(	unsigned int	initialPos,
		const T &	minHorizontalValue,
		const T &	maxHorizontalValue,
		std::vector< T > &	gridValues,
		std::vector< unsigned int > &	bins
	)		const

Calculates the weighted histogram of the sub-sequence.

In this method, gridValues is a std::vector<T>, but internally, it will be copied to an structure of the type UniformOneDGrid<T>.

Definition at line 1831 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

 {
   queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
 
   for (unsigned int j = 0; j < bins.size(); ++j) {
     bins[j] = 0;
   }
 
   double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
   double minCenter = minHorizontalValue - horizontalDelta/2.;
   double maxCenter = maxHorizontalValue + horizontalDelta/2.;
   UniformOneDGrid<T> tmpGrid(m_env,
                                     "",
                                     bins.size(),
                               minCenter,
                               maxCenter);
   gridValues.clear();
   gridValues.resize(tmpGrid.size(),0.);
   for (unsigned int i = 0; i < tmpGrid.size(); ++i) {
     gridValues[i] = tmpGrid[i];
   }
 
   unsigned int dataSize = this->subSequenceSize();
   for (unsigned int j = 0; j < dataSize; ++j) {
     double value = m_seq[j];
     if (value < minHorizontalValue) {
       bins[0] += value;
     }
     else if (value >= maxHorizontalValue) {
       bins[bins.size()-1] += value;
     }
     else {
       unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
       bins[index] += value;
     }
   }
 
   //std::cout << "In ScalarSequence<T>::subWeightHistogram():" << std::endl;
   //for (unsigned int j = 0; j < bins.size(); ++j) {
   //  std::cout << "bins[" << j << "] = " << bins[j] << std::endl;
   //}
 
   return;
 }

template<class T >

void QUESO::ScalarSequence< T >::subWriteContents	(	unsigned int	initialPos,
		unsigned int	numPos,
		const std::string &	fileName,
		const std::string &	fileType,
		const std::set< unsigned int > &	allowedSubEnvIds
	)		const

Writes the sub-sequence to a file.

Given the allowed sub environments (allowedSubEnvIds) that are allowed to write to file, together with the file name and type (fileName, fileType), it writes the entire sub- sequence to the file. The sum of the initial position of the sequence (initialPos) with the number of positions that will be written (numPos) must equal the size of the sequence.

Definition at line 2578 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

Referenced by QUESO::MetropolisHastingsSG< P_V, P_M >::generateFullChain(), and QUESO::MetropolisHastingsSG< P_V, P_M >::generateSequence().

 {
   queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");
 
   FilePtrSetStruct filePtrSet;
   if (m_env.openOutputFile(fileName,
                            fileType,
                            allowedSubEnvIds,
                            false, // A 'true' causes problems when the user chooses (via options
                                   // in the input file) to use just one file for all outputs.
                            filePtrSet)) {
     this->subWriteContents(initialPos,
                            numPos,
                            *filePtrSet.ofsVar,
                            fileType);
     m_env.closeFile(filePtrSet,fileType);
   }
 
   return;
 }

template<class T >

void QUESO::ScalarSequence< T >::subWriteContents	(	unsigned int	initialPos,
		unsigned int	numPos,
		std::ofstream &	ofs,
		const std::string &	fileType
	)		const

Writes the sub-sequence to a file.

Uses object of the type std::ofstream.

Definition at line 2606 of file ScalarSequence.C.

 {
   if (initialPos) {}; // just to remove compiler warning
   if (numPos)     {}; // just to remove compiler warning
   if (&fileType)  {}; // just to remove compiler warning
   ofs << m_name << "_sub" << m_env.subIdString() << " = zeros(" << this->subSequenceSize()
       << ","                                                    << 1
       << ");"
       << std::endl;
   ofs << m_name << "_sub" << m_env.subIdString() << " = [";
   unsigned int chainSize = this->subSequenceSize();
   for (unsigned int j = 0; j < chainSize; ++j) {
     ofs << m_seq[j]
         << std::endl;
   }
   ofs << "];\n";
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::unifiedGaussian1dKde	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		double	unifiedScaleValue,
		const std::vector< T > &	unifiedEvaluationPositions,
		std::vector< double > &	unifiedDensityValues
	)		const

Gaussian kernel for the KDE estimate of the unified sequence.

Definition at line 2349 of file ScalarSequence.C.

References k, QUESO::MiscGaussianDensity(), queso_error_msg, queso_require_msg, RawValue_MPI_DOUBLE, RawValue_MPI_SUM, and RawValue_MPI_UNSIGNED.

Referenced by QUESO::SequenceOfVectors< V, M >::unifiedGaussian1dKde().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subGaussian1dKde(initialPos,
                                   unifiedScaleValue,
                                   unifiedEvaluationPositions,
                                   unifiedDensityValues);
   }
 
   // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
 
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       bool bRC = ((initialPos                        <  this->subSequenceSize()          ) &&
                   (0                                 <  unifiedEvaluationPositions.size()) &&
                   (unifiedEvaluationPositions.size() == unifiedDensityValues.size()      ));
       queso_require_msg(bRC, "invalid input data");
 
       unsigned int localDataSize = this->subSequenceSize() - initialPos;
       unsigned int unifiedDataSize = 0;
       m_env.inter0Comm().Allreduce((void *) &localDataSize, (void *) &unifiedDataSize, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedGaussian1dKde()",
                                    "failed MPI.Allreduce() for data size");
 
       unsigned int numEvals = unifiedEvaluationPositions.size();
 
       std::vector<double> densityValues(numEvals,0.);
       double unifiedScaleInv = 1./unifiedScaleValue;
       for (unsigned int j = 0; j < numEvals; ++j) {
         double x = unifiedEvaluationPositions[j];
         double value = 0.;
         for (unsigned int k = 0; k < localDataSize; ++k) {
           double xk = m_seq[initialPos+k];
           value += MiscGaussianDensity((x-xk)*unifiedScaleInv,0.,1.);
         }
         densityValues[j] = value;
       }
 
       for (unsigned int j = 0; j < numEvals; ++j) {
         unifiedDensityValues[j] = 0.;
       }
       m_env.inter0Comm().Allreduce((void *) &densityValues[0], (void *) &unifiedDensityValues[0], (int) numEvals, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedGaussian1dKde()",
                                    "failed MPI.Allreduce() for density values");
 
       for (unsigned int j = 0; j < numEvals; ++j) {
         unifiedDensityValues[j] *= unifiedScaleInv/((double) unifiedDataSize);
       }
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedGaussian1dKde()"
                                 << ": unifiedDensityValues[0] = "                                       << unifiedDensityValues[0]
                                 << ", unifiedDensityValues[" << unifiedDensityValues.size()-1 << "] = " << unifiedDensityValues[unifiedDensityValues.size()-1]
                                 << std::endl;
       }
     }
     else {
       // Node not in the 'inter0' communicator
       this->subGaussian1dKde(initialPos,
                              unifiedScaleValue,
                              unifiedEvaluationPositions,
                              unifiedDensityValues);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::unifiedHistogram	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		const T &	unifiedMinHorizontalValue,
		const T &	unifiedMaxHorizontalValue,
		std::vector< T > &	unifiedCenters,
		std::vector< unsigned int > &	unifiedBins
	)		const

Calculates the histogram of the unified sequence.

Definition at line 1657 of file ScalarSequence.C.

References queso_error_msg, queso_require_equal_to_msg, queso_require_greater_equal_msg, RawValue_MPI_SUM, and RawValue_MPI_UNSIGNED.

Referenced by QUESO::SequenceOfVectors< V, M >::unifiedHistogram().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subHistogram(initialPos,
                               unifiedMinHorizontalValue,
                               unifiedMaxHorizontalValue,
                               unifiedCenters,
                               unifiedBins);
   }
 
   // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
 
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       queso_require_equal_to_msg(unifiedCenters.size(), unifiedBins.size(), "vectors 'unifiedCenters' and 'unifiedBins' have different sizes");
 
       queso_require_greater_equal_msg(unifiedBins.size(), 3, "number of 'unifiedBins' is too small: should be at least 3");
 
       for (unsigned int j = 0; j < unifiedBins.size(); ++j) {
         unifiedCenters[j] = 0.;
         unifiedBins[j] = 0;
       }
 
       double unifiedHorizontalDelta = (unifiedMaxHorizontalValue - unifiedMinHorizontalValue)/(((double) unifiedBins.size()) - 2.); // IMPORTANT: -2
 
       double unifiedMinCenter = unifiedMinHorizontalValue - unifiedHorizontalDelta/2.;
       double unifiedMaxCenter = unifiedMaxHorizontalValue + unifiedHorizontalDelta/2.;
       for (unsigned int j = 0; j < unifiedCenters.size(); ++j) {
         double factor = ((double) j)/(((double) unifiedCenters.size()) - 1.);
         unifiedCenters[j] = (1. - factor) * unifiedMinCenter + factor * unifiedMaxCenter;
       }
 
       std::vector<unsigned int> localBins(unifiedBins.size(),0);
       unsigned int dataSize = this->subSequenceSize();
       for (unsigned int j = 0; j < dataSize; ++j) {
         double value = m_seq[j];
         if (value < unifiedMinHorizontalValue) {
           localBins[0]++;
         }
         else if (value >= unifiedMaxHorizontalValue) {
           localBins[localBins.size()-1]++;
         }
         else {
           unsigned int index = 1 + (unsigned int) ((value - unifiedMinHorizontalValue)/unifiedHorizontalDelta);
           localBins[index]++;
         }
       }
 
       m_env.inter0Comm().Allreduce((void *) &localBins[0], (void *) &unifiedBins[0], (int) localBins.size(), RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedHistogram()",
                                    "failed MPI.Allreduce() for bins");
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
         for (unsigned int i = 0; i < unifiedCenters.size(); ++i) {
           *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedHistogram()"
                                   << ": i = " << i
                                   << ", unifiedMinHorizontalValue = " << unifiedMinHorizontalValue
                                   << ", unifiedMaxHorizontalValue = " << unifiedMaxHorizontalValue
                                   << ", unifiedCenters = "            << unifiedCenters[i]
                                   << ", unifiedBins = "               << unifiedBins[i]
                                   << std::endl;
         }
       }
     }
     else {
       // Node not in the 'inter0' communicator
       this->subHistogram(initialPos,
                          unifiedMinHorizontalValue,
                          unifiedMaxHorizontalValue,
                          unifiedCenters,
                          unifiedBins);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return;
 }

template<class T >

T QUESO::ScalarSequence< T >::unifiedInterQuantileRange	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos
	)		const

Returns the interquartile range of the values in the unified sequence.

Definition at line 2117 of file ScalarSequence.C.

References queso_error_msg, queso_require_equal_to_msg, RawValue_MPI_SUM, RawValue_MPI_UNSIGNED, and QUESO::ScalarSequence< T >::subSequenceSize().

Referenced by QUESO::SequenceOfVectors< V, M >::unifiedInterQuantileRange().

 {
   T unifiedIqrValue = 0.;
 
   if (m_env.numSubEnvironments() == 1) {
     return this->subInterQuantileRange(initialPos);
   }
 
   // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
 
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedInterQuantileRange(), beginning logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
 
       ScalarSequence unifiedSortedSequence(m_env,0,"");
       this->unifiedSort(useOnlyInter0Comm,
                         initialPos,
                         unifiedSortedSequence);
       unsigned int unifiedDataSize = unifiedSortedSequence.subSequenceSize();
 
       unsigned int localDataSize = this->subSequenceSize() - initialPos;
       unsigned int sumOfLocalSizes = 0;
       m_env.inter0Comm().Allreduce((void *) &localDataSize, (void *) &sumOfLocalSizes, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedInterQuantileRange()",
                                    "failed MPI.Allreduce() for data size");
 
       queso_require_equal_to_msg(sumOfLocalSizes, unifiedDataSize, "incompatible unified sizes");
 
       unsigned int pos1 = (unsigned int) ( (((double) unifiedDataSize) + 1.)*1./4. - 1. );
       unsigned int pos3 = (unsigned int) ( (((double) unifiedDataSize) + 1.)*3./4. - 1. );
 
       double fraction1 = (((double) unifiedDataSize) + 1.)*1./4. - 1. - ((double) pos1);
       double fraction3 = (((double) unifiedDataSize) + 1.)*3./4. - 1. - ((double) pos3);
 
       T value1 = (1.-fraction1) * unifiedSortedSequence[pos1] + fraction1 * unifiedSortedSequence[pos1+1];
       T value3 = (1.-fraction3) * unifiedSortedSequence[pos3] + fraction3 * unifiedSortedSequence[pos3+1];
       unifiedIqrValue = value3 - value1;
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedInterQuantileRange()"
                                 << ": unifiedIqrValue = " << unifiedIqrValue
                                 << ", localDataSize = "   << localDataSize
                                 << ", unifiedDataSize = " << unifiedDataSize
                                 << ", pos1 = "            << pos1
                                 << ", pos3 = "            << pos3
                                 << ", value1 = "          << value1
                                 << ", value3 = "          << value3
                                 << std::endl;
 
         // Save data only once into a separate file
   //std::ofstream* ofsvar = new std::ofstream(("unif_sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
         //if ((ofsvar            == NULL ) ||
         //    (ofsvar->is_open() == false)) {
         //  delete ofsvar;
         //  ofsvar = new std::ofstream(("unif_sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
 
         //  *ofsvar << "var_unif_sort_sub" << m_env.subIdString() << " = zeros(" << 1
         //          << ","                                                       << unifiedDataSize
         //          << ");"
         //          << std::endl;
         //  for (unsigned int j = 0; j < unifiedDataSize; ++j) {
         //    *ofsvar << "var_unif_sort_sub" << m_env.subIdString() << "(" << 1
         //            << ","                                               << j+1
         //            << ") = "                                            << unifiedSortedSequence[j]
         //            << ";"
         //            << std::endl;
         //  }
         //}
         //delete ofsvar;
       }
 
       //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedInterQuantileRange(), ending logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
     }
     else {
       // Node not in the 'inter0' communicator
       unifiedIqrValue = this->subInterQuantileRange(initialPos);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return unifiedIqrValue;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::unifiedMaxPlain ( bool useOnlyInter0Comm ) const

Finds the maximum value of the unified sequence of scalars.

Definition at line 380 of file ScalarSequence.C.

Referenced by QUESO::MLSampling< P_V, P_M >::generateSequence_Step03_inter0().

 {
   if (m_unifiedMaxPlain == NULL) {
     if (m_unifiedMinPlain == NULL) m_unifiedMinPlain = new T(0.);
     m_unifiedMaxPlain = new T(0.);
     unifiedMinMaxExtra(useOnlyInter0Comm,0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
   }
 
   return *m_unifiedMaxPlain;
 }

template<class T >

T QUESO::ScalarSequence< T >::unifiedMeanExtra	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		unsigned int	localNumPos
	)		const

Finds the mean value of the unified sequence of numPos positions starting at position initialPos.

Definition at line 874 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, RawValue_MPI_DOUBLE, RawValue_MPI_SUM, and RawValue_MPI_UNSIGNED.

Referenced by QUESO::ComputeCovCorrBetweenScalarSequences(), QUESO::MLSampling< P_V, P_M >::generateSequence(), and QUESO::SequenceOfVectors< V, M >::unifiedMeanExtra().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subMeanExtra(initialPos,
                               numPos);
   }
 
   // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
 
   T unifiedMeanValue = 0.;
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       bool bRC = ((initialPos          <  this->subSequenceSize()) &&
                   (0                   <  numPos                 ) &&
                   ((initialPos+numPos) <= this->subSequenceSize()));
       queso_require_msg(bRC, "invalid input data");
 
       unsigned int finalPosPlus1 = initialPos + numPos;
       T localSum = 0.;
       for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
         localSum += m_seq[j];
       }
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMeanExtra()"
                                 << ": initialPos = " << initialPos
                                 << ", numPos = "     << numPos
                                 << ", before MPI.Allreduce"
                                 << std::endl;
       }
       //std::cout << m_env.inter0Comm().MyPID()
       //          << std::endl;
       //sleep(1);
       unsigned int unifiedNumPos = 0;
       m_env.inter0Comm().Allreduce((void *) &numPos, (void *) &unifiedNumPos, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedMeanExtra()",
                                    "failed MPI.Allreduce() for numPos");
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMeanExtra()"
                                 << ": numPos = "        << numPos
                                 << ", unifiedNumPos = " << unifiedNumPos
                                 << std::endl;
       }
 
       m_env.inter0Comm().Allreduce((void *) &localSum, (void *) &unifiedMeanValue, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedMeanExtra()",
                                    "failed MPI.Allreduce() for sum");
 
       unifiedMeanValue /= ((T) unifiedNumPos);
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMeanExtra()"
                                 << ": localSum = "         << localSum
                                 << ", unifiedMeanValue = " << unifiedMeanValue
                                 << std::endl;
       }
     }
     else {
       // Node not in the 'inter0' communicator
       this->subMeanExtra(initialPos,
                          numPos);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return unifiedMeanValue;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::unifiedMeanPlain ( bool useOnlyInter0Comm ) const

Finds the mean value of the unified sequence of scalars.

Definition at line 405 of file ScalarSequence.C.

Referenced by QUESO::MLSampling< P_V, P_M >::generateSequence().

 {
   if (m_unifiedMeanPlain == NULL) {
     m_unifiedMeanPlain = new T(0.);
     *m_unifiedMeanPlain = unifiedMeanExtra(useOnlyInter0Comm,0,subSequenceSize());
   }
 
   return *m_unifiedMeanPlain;
 }

template<class T >

T QUESO::ScalarSequence< T >::unifiedMedianExtra	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		unsigned int	localNumPos
	)		const

Finds the median value of the unified sequence, considering numPos positions starting at position initialPos.

Definition at line 995 of file ScalarSequence.C.

References queso_error_msg, and queso_require_msg.

Referenced by QUESO::SequenceOfVectors< V, M >::unifiedMedianExtra().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subMedianExtra(initialPos,
                                 numPos);
   }
 
   // As of 07/Jul/2012, this routine does *not* require sub sequences to have equal size. Good.
 
   T unifiedMedianValue = 0.;
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       bool bRC = ((initialPos          <  this->subSequenceSize()) &&
                   (0                   <  numPos                 ) &&
                   ((initialPos+numPos) <= this->subSequenceSize()));
       queso_require_msg(bRC, "invalid input data");
 
       ScalarSequence unifiedSortedSequence(m_env,0,"");
       this->unifiedSort(useOnlyInter0Comm,
                         initialPos,
                         unifiedSortedSequence);
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMedianExtra()"
                                 << ", unifiedMedianValue = " << unifiedMedianValue
                                 << std::endl;
       }
     }
     else {
       // Node not in the 'inter0' communicator
       this->subMedianExtra(initialPos,
                            numPos);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return unifiedMedianValue;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::unifiedMedianPlain ( bool useOnlyInter0Comm ) const

Finds the median value of the unified sequence of scalars.

Definition at line 429 of file ScalarSequence.C.

 {
   if (m_unifiedMedianPlain == NULL) {
     m_unifiedMedianPlain = new T(0.);
     *m_unifiedMedianPlain = unifiedMedianExtra(useOnlyInter0Comm,0,subSequenceSize());
   }
 
   return *m_unifiedMedianPlain;
 }

template<class T>

void QUESO::ScalarSequence< T >::unifiedMinMaxExtra	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		unsigned int	numPos,
		T &	unifiedMinValue,
		T &	unifiedMaxValue
	)		const

Finds the minimum and the maximum values of the unified sequence, considering numPos positions starting at position initialPos.

Definition at line 1539 of file ScalarSequence.C.

References queso_error_msg, RawValue_MPI_DOUBLE, RawValue_MPI_MAX, and RawValue_MPI_MIN.

Referenced by QUESO::MLSampling< P_V, P_M >::generateSequence(), and QUESO::SequenceOfVectors< V, M >::unifiedMinMaxExtra().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subMinMaxExtra(initialPos,
                                 numPos,
                                 unifiedMinValue,
                                 unifiedMaxValue);
   }
 
   // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
 
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       // Find local min and max
       T minValue;
       T maxValue;
       this->subMinMaxExtra(initialPos,
                            numPos,
                            minValue,
                            maxValue);
 
       // Get overall min
       std::vector<double> sendBuf(1,0.);
       for (unsigned int i = 0; i < sendBuf.size(); ++i) {
         sendBuf[i] = minValue;
       }
       m_env.inter0Comm().Allreduce((void *) &sendBuf[0], (void *) &unifiedMinValue, (int) sendBuf.size(), RawValue_MPI_DOUBLE, RawValue_MPI_MIN,
                                    "ScalarSequence<T>::unifiedMinMaxExtra()",
                                    "failed MPI.Allreduce() for min");
 
       // Get overall max
       for (unsigned int i = 0; i < sendBuf.size(); ++i) {
         sendBuf[i] = maxValue;
       }
       m_env.inter0Comm().Allreduce((void *) &sendBuf[0], (void *) &unifiedMaxValue, (int) sendBuf.size(), RawValue_MPI_DOUBLE, RawValue_MPI_MAX,
                                    "ScalarSequence<T>::unifiedMinMaxExtra()",
                                    "failed MPI.Allreduce() for max");
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMinMaxExtra()"
                                 << ": localMinValue = "   << minValue
                                 << ", localMaxValue = "   << maxValue
                                 << ", unifiedMinValue = " << unifiedMinValue
                                 << ", unifiedMaxValue = " << unifiedMaxValue
                                 << std::endl;
       }
     }
     else {
       // Node not in the 'inter0' communicator
       this->subMinMaxExtra(initialPos,
                            numPos,
                            unifiedMinValue,
                            unifiedMaxValue);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::unifiedMinPlain ( bool useOnlyInter0Comm ) const

Finds the minimum value of the unified sequence of scalars.

Definition at line 354 of file ScalarSequence.C.

 {
   if (m_unifiedMinPlain == NULL) {
     m_unifiedMinPlain = new T(0.);
     if (m_unifiedMaxPlain == NULL) m_unifiedMaxPlain = new T(0.);
     unifiedMinMaxExtra(useOnlyInter0Comm,0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
   }
 
   return *m_unifiedMinPlain;
 }

template<class T>

T QUESO::ScalarSequence< T >::unifiedPopulationVariance	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		unsigned int	numPos,
		const T &	unifiedMeanValue
	)		const

Finds the population variance of the unified sequence, considering numPos positions starting at position initialPos and of mean meanValue.

Output: calculated population variance of the unified sequence of scalars.

Definition at line 1243 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, RawValue_MPI_DOUBLE, RawValue_MPI_SUM, and RawValue_MPI_UNSIGNED.

Referenced by QUESO::SequenceOfVectors< V, M >::unifiedPopulationVariance().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subPopulationVariance(initialPos,
                                        numPos,
                                        unifiedMeanValue);
   }
 
   // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
 
   T unifiedPopValue = 0.;
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       bool bRC = ((initialPos          <  this->subSequenceSize()) &&
                   (0                   <  numPos                 ) &&
                   ((initialPos+numPos) <= this->subSequenceSize()));
       queso_require_msg(bRC, "invalid input data");
 
       unsigned int finalPosPlus1 = initialPos + numPos;
       T diff;
       T localPopValue = 0.;
       for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
         diff = m_seq[j] - unifiedMeanValue;
         localPopValue += diff*diff;
       }
 
       unsigned int unifiedNumPos = 0;
       m_env.inter0Comm().Allreduce((void *) &numPos, (void *) &unifiedNumPos, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedPopulationVariance()",
                                    "failed MPI.Allreduce() for numPos");
 
       m_env.inter0Comm().Allreduce((void *) &localPopValue, (void *) &unifiedPopValue, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedPopulationVariance()",
                                    "failed MPI.Allreduce() for popValue");
 
       unifiedPopValue /= ((T) unifiedNumPos);
     }
     else {
       // Node not in the 'inter0' communicator
       this->subPopulationVariance(initialPos,
                                   numPos,
                                   unifiedMeanValue);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return unifiedPopValue;
 }

template<class T>

T QUESO::ScalarSequence< T >::unifiedPositionsOfMaximum	(	const ScalarSequence< T > &	subCorrespondingScalarValues,
		ScalarSequence< T > &	unifiedPositionsOfMaximum
	)

Finds the positions where the maximum element occurs in the unified sequence.

Definition at line 2548 of file ScalarSequence.C.

References queso_require_equal_to_msg, QUESO::ScalarSequence< T >::resizeSequence(), QUESO::ScalarSequence< T >::subMaxPlain(), and QUESO::ScalarSequence< T >::subSequenceSize().

 {
   queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
 
   T maxValue = subCorrespondingScalarValues.subMaxPlain();
   unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
 
   unsigned int numPos = 0;
   for (unsigned int i = 0; i < iMax; ++i) {
     if (subCorrespondingScalarValues[i] == maxValue) {
       numPos++;
     }
   }
 
   unifiedPositionsOfMaximum.resizeSequence(numPos);
   unsigned int j = 0;
   for (unsigned int i = 0; i < iMax; ++i) {
     if (subCorrespondingScalarValues[i] == maxValue) {
       unifiedPositionsOfMaximum[j] = (*this)[i];
       j++;
     }
   }
 
   return maxValue;
 }

template<class T >

void QUESO::ScalarSequence< T >::unifiedReadContents	(	const std::string &	fileName,
		const std::string &	fileType,
		const unsigned int	subSequenceSize
	)

Reads the unified sequence from a file.

Definition at line 2833 of file ScalarSequence.C.

References dataIn, QUESO::FilePtrSetStruct::ifsVar, QUESO::MiscGetEllapsedSeconds(), queso_error_msg, queso_require_equal_to_msg, queso_require_greater_equal_msg, queso_require_less_msg, UQ_FILE_EXTENSION_FOR_HDF_FORMAT, UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, and QUESO::UQ_OK_RC.

Referenced by QUESO::MLSampling< P_V, P_M >::restartML().

 {
   std::string fileType(inputFileType);
 #ifdef QUESO_HAS_HDF5
   // Do nothing
 #else
   if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
     if (m_env.subDisplayFile()) {
       *m_env.subDisplayFile() << "WARNING in ScalarSequence<T>::unifiedReadContents()"
                               << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
                               << "' has been requested, but this QUESO library has not been built with 'hdf5'"
                               << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
                               << "' instead..."
                               << std::endl;
     }
     if (m_env.subRank() == 0) {
       std::cerr << "WARNING in ScalarSequence<T>::unifiedReadContents()"
                 << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
                 << "' has been requested, but this QUESO library has not been built with 'hdf5'"
                 << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
                 << "' instead..."
                 << std::endl;
     }
     fileType = UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT;
   }
 #endif
 
   //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
   if (m_env.subDisplayFile()) {
     *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedReadContents()"
                             << ": worldRank "                      << m_env.worldRank()
                             << ", fullRank "                       << m_env.fullRank()
                             << ", subEnvironment "                 << m_env.subId()
                             << ", subRank "                        << m_env.subRank()
                             << ", inter0Rank "                     << m_env.inter0Rank()
       //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
                             << ", fileName = "                     << fileName
                             << ", subReadSize = "                  << subReadSize
       //<< ", unifiedReadSize = "              << unifiedReadSize
                             << std::endl;
   }
 
   this->resizeSequence(subReadSize);
 
   if (m_env.inter0Rank() >= 0) {
     double unifiedReadSize = subReadSize*m_env.inter0Comm().NumProc();
 
     // In the logic below, the id of a line' begins with value 0 (zero)
     unsigned int idOfMyFirstLine = 1 + m_env.inter0Rank()*subReadSize;
     unsigned int idOfMyLastLine = (1 + m_env.inter0Rank())*subReadSize;
     unsigned int numParams = 1; // this->vectorSizeLocal();
 
     for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) { // "m or hdf"
       if (m_env.inter0Rank() == (int) r) {
         // My turn
         FilePtrSetStruct unifiedFilePtrSet;
         if (m_env.openUnifiedInputFile(fileName,
                                        fileType,
                                        unifiedFilePtrSet)) {
           if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
             if (r == 0) {
               // Read number of chain positions in the file by taking care of the first line,
               // which resembles something like 'variable_name = zeros(n_positions,m_params);'
         std::string tmpString;
 
               // Read 'variable name' string
               *unifiedFilePtrSet.ifsVar >> tmpString;
         //std::cout << "Just read '" << tmpString << "'" << std::endl;
 
               // Read '=' sign
               *unifiedFilePtrSet.ifsVar >> tmpString;
           //std::cout << "Just read '" << tmpString << "'" << std::endl;
               queso_require_equal_to_msg(tmpString, "=", "string should be the '=' sign");
 
               // Read 'zeros(n_positions,n_params)' string
               *unifiedFilePtrSet.ifsVar >> tmpString;
         //std::cout << "Just read '" << tmpString << "'" << std::endl;
               unsigned int posInTmpString = 6;
 
               // Isolate 'n_positions' in a string
               //char nPositionsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
         std::string nPositionsString((size_t) (tmpString.size()-posInTmpString+1),' ');
               unsigned int posInPositionsString = 0;
               do {
                 queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ',' not found in first line of file");
                 nPositionsString[posInPositionsString++] = tmpString[posInTmpString++];
               } while (tmpString[posInTmpString] != ',');
               nPositionsString[posInPositionsString] = '\0';
 
               // Isolate 'n_params' in a string
               posInTmpString++; // Avoid reading ',' char
               //char nParamsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
         std::string nParamsString((size_t) (tmpString.size()-posInTmpString+1),' ');
               unsigned int posInParamsString = 0;
               do {
                 queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ')' not found in first line of file");
                 nParamsString[posInParamsString++] = tmpString[posInTmpString++];
               } while (tmpString[posInTmpString] != ')');
               nParamsString[posInParamsString] = '\0';
 
               // Convert 'n_positions' and 'n_params' strings to numbers
               unsigned int sizeOfChainInFile = (unsigned int) strtod(nPositionsString.c_str(),NULL);
               unsigned int numParamsInFile   = (unsigned int) strtod(nParamsString.c_str(),   NULL);
               if (m_env.subDisplayFile()) {
                 *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedReadContents()"
                                         << ": worldRank "           << m_env.worldRank()
                                         << ", fullRank "            << m_env.fullRank()
                                         << ", sizeOfChainInFile = " << sizeOfChainInFile
                                         << ", numParamsInFile = "   << numParamsInFile
                                         << std::endl;
               }
 
               // Check if [size of chain in file] >= [requested unified sequence size]
               queso_require_greater_equal_msg(sizeOfChainInFile, unifiedReadSize, "size of chain in file is not big enough");
 
               // Check if [num params in file] == [num params in current chain]
               queso_require_equal_to_msg(numParamsInFile, numParams, "number of parameters of chain in file is different than number of parameters in this chain object");
             } // if (r == 0)
 
             // Code common to any core in 'inter0Comm', including core of rank 0
             unsigned int maxCharsPerLine = 64*numParams; // Up to about 60 characters to represent each parameter value
 
             unsigned int lineId = 0;
             while (lineId < idOfMyFirstLine) {
               unifiedFilePtrSet.ifsVar->ignore(maxCharsPerLine,'\n');
               lineId++;
             };
 
             if (r == 0) {
               // Take care of initial part of the first data line,
               // which resembles something like 'variable_name = [value1 value2 ...'
         std::string tmpString;
 
               // Read 'variable name' string
               *unifiedFilePtrSet.ifsVar >> tmpString;
           //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
 
               // Read '=' sign
               *unifiedFilePtrSet.ifsVar >> tmpString;
         //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
               queso_require_equal_to_msg(tmpString, "=", "in core 0, string should be the '=' sign");
 
               // Take into account the ' [' portion
         std::streampos tmpPos = unifiedFilePtrSet.ifsVar->tellg();
               unifiedFilePtrSet.ifsVar->seekg(tmpPos+(std::streampos)2);
             }
 
             T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
             while (lineId <= idOfMyLastLine) {
               for (unsigned int i = 0; i < numParams; ++i) {
                 *unifiedFilePtrSet.ifsVar >> tmpScalar;
               }
               m_seq[lineId - idOfMyFirstLine] = tmpScalar;
               lineId++;
             };
           }
 #ifdef QUESO_HAS_HDF5
           else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
             if (r == 0) {
               hid_t dataset = H5Dopen2(unifiedFilePtrSet.h5Var,
                                        "seq_of_vectors",
                                        H5P_DEFAULT); // Dataset access property list
               hid_t datatype  = H5Dget_type(dataset);
               H5T_class_t t_class = H5Tget_class(datatype);
               queso_require_equal_to_msg(t_class, H5T_FLOAT, "t_class is not H5T_DOUBLE");
               hid_t dataspace = H5Dget_space(dataset);
               int   rank      = H5Sget_simple_extent_ndims(dataspace);
               queso_require_equal_to_msg(rank, 2, "hdf rank is not 2");
               hsize_t dims_in[2];
               int     status_n;
               status_n  = H5Sget_simple_extent_dims(dataspace, dims_in, NULL);
               if (status_n) {}; // just to remove compiler warning
         //std::cout << "In ScalarSequence<T>::unifiedReadContents()"
               //          << ": dims_in[0] = " << dims_in[0]
               //          << ", dims_in[1] = " << dims_in[1]
               //          << std::endl;
               queso_require_equal_to_msg(dims_in[0], numParams, "dims_in[0] is not equal to 'numParams'");
               queso_require_greater_equal_msg(dims_in[1], subReadSize, "dims_in[1] is smaller that requested 'subReadSize'");
 
               struct timeval timevalBegin;
               int iRC = UQ_OK_RC;
               iRC = gettimeofday(&timevalBegin,NULL);
               if (iRC) {}; // just to remove compiler warning
 
               unsigned int chainSizeIn = (unsigned int) dims_in[1];
               //double* dataIn[numParams]; // avoid compiler warning
         std::vector<double*> dataIn((size_t) numParams,NULL);
               dataIn[0] = (double*) malloc(numParams*chainSizeIn*sizeof(double));
               for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1'
                 dataIn[i] = dataIn[i-1] + chainSizeIn; // Yes, just 'chainSizeIn', not 'chainSizeIn*sizeof(double)'
               }
               //std::cout << "In ScalarSequence<T>::unifiedReadContents(): h5 case, memory allocated" << std::endl;
               herr_t status;
               status = H5Dread(dataset,
                                H5T_NATIVE_DOUBLE,
                                H5S_ALL,
                                dataspace,
                                H5P_DEFAULT,
                                dataIn[0]);
               if (status) {}; // just to remove compiler warning
               //std::cout << "In ScalarSequence<T>::unifiedReadContents(): h5 case, data read" << std::endl;
               T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
               for (unsigned int j = 0; j < subReadSize; ++j) { // Yes, 'subReadSize', not 'chainSizeIn'
                 for (unsigned int i = 0; i < numParams; ++i) {
                   tmpScalar = dataIn[i][j];
                 }
                 m_seq[j] = tmpScalar;
               }
 
               double readTime = MiscGetEllapsedSeconds(&timevalBegin);
               if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
                 *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedReadContents()"
                                         << ": worldRank "      << m_env.worldRank()
                                         << ", fullRank "       << m_env.fullRank()
                                         << ", subEnvironment " << m_env.subId()
                                         << ", subRank "        << m_env.subRank()
                                         << ", inter0Rank "     << m_env.inter0Rank()
                                         << ", fileName = "     << fileName
                                         << ", numParams = "    << numParams
                                         << ", chainSizeIn = "  << chainSizeIn
                                         << ", subReadSize = "  << subReadSize
                                         << ", readTime = "     << readTime << " seconds"
                                         << std::endl;
               }
 
               H5Sclose(dataspace);
               H5Tclose(datatype);
               H5Dclose(dataset);
               //free(dataIn[0]); // related to the changes above for compiler warning
               for (unsigned int tmpIndex = 0; tmpIndex < dataIn.size(); tmpIndex++) {
                 free (dataIn[tmpIndex]);
               }
             }
             else {
               queso_error_msg("hdf file type not supported for multiple sub-environments yet");
             }
           }
 #endif
           else {
             queso_error_msg("invalid file type");
           }
           m_env.closeFile(unifiedFilePtrSet,fileType);
         } // if (m_env.openUnifiedInputFile())
       } // if (m_env.inter0Rank() == (int) r)
       m_env.inter0Comm().Barrier();
     } // for r
   } // if (m_env.inter0Rank() >= 0)
   else {
     T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
     for (unsigned int i = 1; i < subReadSize; ++i) {
       m_seq[i] = tmpScalar;
     }
   }
 
   if (m_env.subDisplayFile()) {
     *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedReadContents()"
                             << ", fileName = " << fileName
                             << std::endl;
   }
   //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
 
   return;
 }

template<class T>

T QUESO::ScalarSequence< T >::unifiedSampleStd	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		unsigned int	localNumPos,
		const T &	unifiedMeanValue
	)		const

Finds the sample standard deviation of the unified sequence, considering localnumPos positions starting at position initialPos and of mean unifiedMeanValue.

Definition at line 1156 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, RawValue_MPI_DOUBLE, RawValue_MPI_SUM, and RawValue_MPI_UNSIGNED.

Referenced by QUESO::SequenceOfVectors< V, M >::unifiedSampleStd().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subSampleStd(initialPos,
                                    numPos,
                                    unifiedMeanValue);
   }
 
   // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
 
   T unifiedStdValue = 0.;
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       bool bRC = ((initialPos          <  this->subSequenceSize()) &&
                   (0                   <  numPos                 ) &&
                   ((initialPos+numPos) <= this->subSequenceSize()));
       queso_require_msg(bRC, "invalid input data");
 
       unsigned int finalPosPlus1 = initialPos + numPos;
       T diff;
       T localStdValue = 0.;
       for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
         diff = m_seq[j] - unifiedMeanValue;
         localStdValue += diff*diff;
       }
 
       unsigned int unifiedNumPos = 0;
       m_env.inter0Comm().Allreduce((void *) &numPos, (void *) &unifiedNumPos, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedSampleStd()",
                                    "failed MPI.Allreduce() for numPos");
 
       m_env.inter0Comm().Allreduce((void *) &localStdValue, (void *) &unifiedStdValue, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedSampleStd()",
                                    "failed MPI.Allreduce() for stdValue");
 
       unifiedStdValue /= (((T) unifiedNumPos) - 1.);
       unifiedStdValue = sqrt(unifiedStdValue);
     }
     else {
       // Node not in the 'inter0' communicator
       this->subSampleStd(initialPos,
                          numPos,
                          unifiedMeanValue);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return unifiedStdValue;
 }

template<class T>

T QUESO::ScalarSequence< T >::unifiedSampleVarianceExtra	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		unsigned int	localNumPos,
		const T &	unifiedMeanValue
	)		const

Finds the sample variance of the unified sequence, considering numPos positions starting at position initialPos and of mean meanVec.

Output:

Parameters

samVec is the vector of the calculated sample variance of the unified sequence of scalars.

Definition at line 1069 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, RawValue_MPI_DOUBLE, RawValue_MPI_SUM, and RawValue_MPI_UNSIGNED.

Referenced by QUESO::ComputeCovCorrBetweenScalarSequences(), and QUESO::SequenceOfVectors< V, M >::unifiedSampleVarianceExtra().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subSampleVarianceExtra(initialPos,
                                         numPos,
                                         unifiedMeanValue);
   }
 
   // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
 
   T unifiedSamValue = 0.;
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       bool bRC = ((initialPos          <  this->subSequenceSize()) &&
                   (0                   <  numPos                 ) &&
                   ((initialPos+numPos) <= this->subSequenceSize()));
       queso_require_msg(bRC, "invalid input data");
 
       unsigned int finalPosPlus1 = initialPos + numPos;
       T diff;
       T localSamValue = 0.;
       for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
         diff = m_seq[j] - unifiedMeanValue;
         localSamValue += diff*diff;
       }
 
       unsigned int unifiedNumPos = 0;
       m_env.inter0Comm().Allreduce((void *) &numPos, (void *) &unifiedNumPos, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedSampleVarianceExtra()",
                                    "failed MPI.Allreduce() for numPos");
 
       m_env.inter0Comm().Allreduce((void *) &localSamValue, (void *) &unifiedSamValue, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedSampleVarianceExtra()",
                                    "failed MPI.Allreduce() for samValue");
 
       unifiedSamValue /= (((T) unifiedNumPos) - 1.);
     }
     else {
       // Node not in the 'inter0' communicator
       this->subSampleVarianceExtra(initialPos,
                                    numPos,
                                    unifiedMeanValue);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return unifiedSamValue;
 }

template<class T >

const T & QUESO::ScalarSequence< T >::unifiedSampleVariancePlain ( bool useOnlyInter0Comm ) const

Finds the variance of a sample of the unified sequence of scalars.

Definition at line 453 of file ScalarSequence.C.

 {
   if (m_unifiedSampleVariancePlain == NULL) {
     m_unifiedSampleVariancePlain = new T(0.);
     *m_unifiedSampleVariancePlain = unifiedSampleVarianceExtra(useOnlyInter0Comm,0,subSequenceSize(),unifiedMeanPlain(useOnlyInter0Comm));
   }
 
   return *m_unifiedSampleVariancePlain;
 }

template<class T>

T QUESO::ScalarSequence< T >::unifiedScaleForKde	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		const T &	unifiedIqrValue,
		unsigned int	kdeDimension
	)		const

Selects the scales (bandwidth) for the kernel density estimation, considering the unified sequence.

Definition at line 2248 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, RawValue_MPI_SUM, and RawValue_MPI_UNSIGNED.

Referenced by QUESO::SequenceOfVectors< V, M >::unifiedScalesForKde().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subScaleForKde(initialPos,
                                 unifiedIqrValue,
                                 kdeDimension);
   }
 
   // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
 
   T unifiedScaleValue = 0.;
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       bool bRC = (initialPos <  this->subSequenceSize());
       queso_require_msg(bRC, "invalid input data");
 
       unsigned int localDataSize = this->subSequenceSize() - initialPos;
 
       T unifiedMeanValue = this->unifiedMeanExtra(useOnlyInter0Comm,
                                                   initialPos,
                                                   localDataSize);
 
       T unifiedSamValue = this->unifiedSampleVarianceExtra(useOnlyInter0Comm,
                                                            initialPos,
                                                            localDataSize,
                                                            unifiedMeanValue);
 
       unsigned int unifiedDataSize = 0;
       m_env.inter0Comm().Allreduce((void *) &localDataSize, (void *) &unifiedDataSize, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedScaleForKde()",
                                    "failed MPI.Allreduce() for data size");
 
       if (unifiedIqrValue <= 0.) {
         unifiedScaleValue = 1.06*std::sqrt(unifiedSamValue)/std::pow(unifiedDataSize,1./(4. + ((double) kdeDimension)));
       }
       else {
         unifiedScaleValue = 1.06*std::min(std::sqrt(unifiedSamValue),unifiedIqrValue/1.34)/std::pow(unifiedDataSize,1./(4. + ((double) kdeDimension)));
       }
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedScaleForKde()"
                                 << ": unifiedIqrValue = "   << unifiedIqrValue
                                 << ", unifiedMeanValue = "  << unifiedMeanValue
                                 << ", unifiedSamValue = "   << unifiedSamValue
                                 << ", unifiedDataSize = "   << unifiedDataSize
                                 << ", unifiedScaleValue = " << unifiedScaleValue
                                 << std::endl;
       }
     }
     else {
       // Node not in the 'inter0' communicator
       unifiedScaleValue = this->subScaleForKde(initialPos,
                                                unifiedIqrValue,
                                                kdeDimension);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return unifiedScaleValue;
 }

template<class T >

unsigned int QUESO::ScalarSequence< T >::unifiedSequenceSize ( bool useOnlyInter0Comm ) const

Size of the unified sequence of scalars.

Definition at line 143 of file ScalarSequence.C.

References queso_error_msg, RawValue_MPI_SUM, and RawValue_MPI_UNSIGNED.

Referenced by QUESO::MLSampling< P_V, P_M >::checkpointML(), QUESO::MLSampling< P_V, P_M >::generateSequence_Step02_inter0(), QUESO::MLSampling< P_V, P_M >::generateSequence_Step03_inter0(), QUESO::MLSampling< P_V, P_M >::generateSequence_Step05_inter0(), and QUESO::MLSampling< P_V, P_M >::generateSequence_Step09_all().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subSequenceSize();
   }
 
   // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
 
   unsigned int unifiedNumSamples = 0;
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       unsigned int subNumSamples = this->subSequenceSize();
       m_env.inter0Comm().Allreduce((void *) &subNumSamples, (void *) &unifiedNumSamples, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedSequenceSize()",
                                    "failed MPI.Allreduce() for unifiedSequenceSize()");
     }
     else {
       // Node not in the 'inter0' communicator
       unifiedNumSamples = this->subSequenceSize();
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   return unifiedNumSamples;
 }

template<class T>

void QUESO::ScalarSequence< T >::unifiedSort	(	bool	useOnlyInter0Comm,
		unsigned int	initialPos,
		ScalarSequence< T > &	unifiedSortedSequence
	)		const

Sorts the unified sequence of scalars.

Definition at line 1900 of file ScalarSequence.C.

References queso_error_msg, queso_require_equal_to_msg, QUESO::ScalarSequence< T >::rawData(), RawValue_MPI_ANY_SOURCE, RawValue_MPI_DOUBLE, RawValue_MPI_SUM, RawValue_MPI_UNSIGNED, QUESO::ScalarSequence< T >::resizeSequence(), SCALAR_SEQUENCE_INIT_MPI_MSG, and QUESO::ScalarSequence< T >::subSequenceSize().

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subSort(initialPos,unifiedSortedSequence);
   }
 
   // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
 
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), beginning logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
 
       unsigned int localNumPos = this->subSequenceSize() - initialPos;
 
       std::vector<T> leafData(localNumPos,0.);
       this->extractRawData(0,
                            1,
                            localNumPos,
                            leafData);
 
       if (m_env.inter0Rank() == 0) {
         int minus1NumTreeLevels = 0;
         int power2NumTreeNodes  = 1;
 
         while (power2NumTreeNodes < m_env.inter0Comm().NumProc()) {
           power2NumTreeNodes += power2NumTreeNodes;
           minus1NumTreeLevels++;
         }
 
         if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
           *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedSort()"
                                   << ": sorting tree has " << m_env.inter0Comm().NumProc()
                                   << " nodes and "         << minus1NumTreeLevels+1
                                   << " levels"
                                   << std::endl;
         }
 
         this->parallelMerge(unifiedSortedSequence.rawData(),
                             leafData,
                             minus1NumTreeLevels);
       }
       else if (m_env.inter0Rank() > 0) { // KAUST
         unsigned int uintBuffer[1];
         RawType_MPI_Status status;
         m_env.inter0Comm().Recv((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, RawValue_MPI_ANY_SOURCE, SCALAR_SEQUENCE_INIT_MPI_MSG, &status,
                                 "ScalarSequence<T>::unifiedSort()",
                                 "failed MPI.Recv() for init");
   //if (status) {}; // just to remove compiler warning
 
         unsigned int treeLevel = uintBuffer[0];
         this->parallelMerge(unifiedSortedSequence.rawData(),
                             leafData,
                             treeLevel);
       }
 
       //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), returned from parallelMerge()",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
 
       // Broadcast
       unsigned int unifiedDataSize = unifiedSortedSequence.subSequenceSize();
       m_env.inter0Comm().Bcast((void *) &unifiedDataSize, (int) 1, RawValue_MPI_UNSIGNED, 0,
                                "ScalarSequence<T>::unifiedSort()",
                                "failed MPI.Bcast() for unified data size");
 
       unsigned int sumOfNumPos = 0;
       m_env.inter0Comm().Allreduce((void *) &localNumPos, (void *) &sumOfNumPos, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
                                    "ScalarSequence<T>::unifiedSort()",
                                    "failed MPI.Allreduce() for data size");
 
       queso_require_equal_to_msg(sumOfNumPos, unifiedDataSize, "incompatible unified sizes");
 
       unifiedSortedSequence.resizeSequence(unifiedDataSize);
       m_env.inter0Comm().Bcast((void *) &unifiedSortedSequence.rawData()[0], (int) unifiedDataSize, RawValue_MPI_DOUBLE, 0,
                                "ScalarSequence<T>::unifiedSort()",
                                "failed MPI.Bcast() for unified data");
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
         *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
                                 << ": tree node "                                                                    << m_env.inter0Rank()
                                 << ", unifiedSortedSequence[0] = "                                                   << unifiedSortedSequence[0]
                                 << ", unifiedSortedSequence[" << unifiedSortedSequence.subSequenceSize()-1 << "] = " << unifiedSortedSequence[unifiedSortedSequence.subSequenceSize()-1]
                                 << std::endl;
       }
 
       //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), ending logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
     }
     else {
       // Node not in the 'inter0' communicator
       this->subSort(initialPos,unifiedSortedSequence);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return;
 }

template<class T>

void QUESO::ScalarSequence< T >::unifiedUniformlySampledCdf	(	bool	useOnlyInter0Comm,
		unsigned int	numIntervals,
		T &	unifiedMinDomainValue,
		T &	unifiedMaxDomainValue,
		std::vector< T > &	unifiedCdfValues
	)		const

Uniformly samples from the CDF from the unified sequence.

Definition at line 611 of file ScalarSequence.C.

References queso_error_msg.

 {
   if (m_env.numSubEnvironments() == 1) {
     return this->subUniformlySampledCdf(numEvaluationPoints,
                                         unifiedMinDomainValue,
                                         unifiedMaxDomainValue,
                                         unifiedCdfValues);
   }
 
   // KAUST2
   // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
 
   if (useOnlyInter0Comm) {
     if (m_env.inter0Rank() >= 0) {
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
         *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedUniformlySampledCdf()"
                                 << std::endl;
       }
 
       T                         unifiedTmpMinValue;
       T                         unifiedTmpMaxValue;
       std::vector<T>            unifiedCenters(numEvaluationPoints,0.);
       std::vector<unsigned int> unifiedBins   (numEvaluationPoints,0);
 
       this->unifiedMinMaxExtra(useOnlyInter0Comm,
                                0, // initialPos
                                this->subSequenceSize(),
                                unifiedTmpMinValue,
                                unifiedTmpMaxValue);
       this->unifiedHistogram(useOnlyInter0Comm,
                              0, // initialPos
                              unifiedTmpMinValue,
                              unifiedTmpMaxValue,
                              unifiedCenters,
                              unifiedBins);
 
       unifiedMinDomainValue = unifiedCenters[0];
       unifiedMaxDomainValue = unifiedCenters[unifiedCenters.size()-1];
 
       unsigned int unifiedTotalSumOfBins = 0;
       for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
         unifiedTotalSumOfBins += unifiedBins[i];
       }
 
       std::vector<unsigned int> unifiedPartialSumsOfBins(numEvaluationPoints,0);
       unifiedPartialSumsOfBins[0] = unifiedBins[0];
       for (unsigned int i = 1; i < numEvaluationPoints; ++i) { // Yes, from '1'
         unifiedPartialSumsOfBins[i] = unifiedPartialSumsOfBins[i-1] + unifiedBins[i];
       }
 
       unifiedCdfValues.clear();
       unifiedCdfValues.resize(numEvaluationPoints);
       for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
         unifiedCdfValues[i] = ((T) unifiedPartialSumsOfBins[i])/((T) unifiedTotalSumOfBins);
       }
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
         for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
           *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedUniformlySampledCdf()"
                                   << ": i = " << i
                                   << ", unifiedTmpMinValue = "       << unifiedTmpMinValue
                                   << ", unifiedTmpMaxValue = "       << unifiedTmpMaxValue
                                   << ", unifiedBins = "              << unifiedBins[i]
                                   << ", unifiedCdfValue = "          << unifiedCdfValues[i]
                                   << ", unifiedPartialSumsOfBins = " << unifiedPartialSumsOfBins[i]
                                   << ", unifiedTotalSumOfBins = "    << unifiedTotalSumOfBins
                                   << std::endl;
         }
       }
 
       if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
         *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedUniformlySampledCdf()"
                                 << std::endl;
       }
     }
     else {
       // Node not in the 'inter0' communicator
       this->subUniformlySampledCdf(numEvaluationPoints,
                                    unifiedMinDomainValue,
                                    unifiedMaxDomainValue,
                                    unifiedCdfValues);
     }
   }
   else {
     queso_error_msg("parallel vectors not supported yet");
   }
 
   //m_env.fullComm().Barrier();
 
   return;
 }

template<class T >

void QUESO::ScalarSequence< T >::unifiedWriteContents	(	const std::string &	fileName,
		const std::string &	fileType
	)		const

Writes the unified sequence to a file.

Writes the unified sequence in Matlab/Octave format or, if enabled, in HDF5 format.

Definition at line 2632 of file ScalarSequence.C.

References QUESO::MiscGetEllapsedSeconds(), QUESO::FilePtrSetStruct::ofsVar, queso_error_msg, UQ_FILE_EXTENSION_FOR_HDF_FORMAT, UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, and QUESO::UQ_OK_RC.

Referenced by QUESO::MLSampling< P_V, P_M >::checkpointML(), QUESO::MLSampling< P_V, P_M >::generateSequence(), QUESO::MetropolisHastingsSG< P_V, P_M >::generateSequence(), QUESO::MLSampling< P_V, P_M >::generateSequence_Level0_all(), and QUESO::MLSampling< P_V, P_M >::generateSequence_Step11_inter0().

 {
   std::string fileType(inputFileType);
 #ifdef QUESO_HAS_HDF5
   // Do nothing
 #else
   if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
     if (m_env.subDisplayFile()) {
       *m_env.subDisplayFile() << "WARNING in ScalarSequence<T>::unifiedWriteContents()"
                               << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
                               << "' has been requested, but this QUESO library has not been built with 'hdf5'"
                               << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
                               << "' instead..."
                               << std::endl;
     }
     if (m_env.subRank() == 0) {
       std::cerr << "WARNING in ScalarSequence<T>::unifiedWriteContents()"
                 << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
                 << "' has been requested, but this QUESO library has not been built with 'hdf5'"
                 << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
                 << "' instead..."
                 << std::endl;
     }
     fileType = UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT;
   }
 #endif
 
   // All processors in 'fullComm' should call this routine...
 
   //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
     *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedWriteContents()"
                             << ": worldRank "      << m_env.worldRank()
                             << ", subEnvironment " << m_env.subId()
                             << ", subRank "        << m_env.subRank()
                             << ", inter0Rank "     << m_env.inter0Rank()
       //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
                             << ", fileName = "     << fileName
                             << ", fileType = "     << fileType
                             << std::endl;
   }
 
   // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
 
   if (m_env.inter0Rank() >= 0) {
     for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) {
       if (m_env.inter0Rank() == (int) r) {
         // My turn
         FilePtrSetStruct unifiedFilePtrSet;
         // bool writeOver = (r == 0);
         bool writeOver = false; // A 'true' causes problems when the user chooses (via options
                                 // in the input file) to use just one file for all outputs.
         //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 000 \n" << std::endl;
         if (m_env.openUnifiedOutputFile(fileName,
                                         fileType, // "m or hdf"
                                         writeOver,
                                         unifiedFilePtrSet)) {
           //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 001 \n" << std::endl;
 
           unsigned int chainSize = this->subSequenceSize();
           if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
             if (r == 0) {
               *unifiedFilePtrSet.ofsVar << m_name << "_unified" << " = zeros(" << this->subSequenceSize()*m_env.inter0Comm().NumProc()
                                         << ","                                 << 1
                                         << ");"
                                         << std::endl;
               *unifiedFilePtrSet.ofsVar << m_name << "_unified" << " = [";
             }
 
             for (unsigned int j = 0; j < chainSize; ++j) {
               *unifiedFilePtrSet.ofsVar << m_seq[j]
                                         << std::endl;
             }
 
             m_env.closeFile(unifiedFilePtrSet,fileType);
     }
 #ifdef QUESO_HAS_HDF5
           else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
             unsigned int numParams = 1; // m_vectorSpace.dimLocal();
             if (r == 0) {
               hid_t datatype = H5Tcopy(H5T_NATIVE_DOUBLE);
               //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data type created" << std::endl;
               hsize_t dimsf[2];
               dimsf[0] = numParams;
               dimsf[1] = chainSize;
               hid_t dataspace = H5Screate_simple(2, dimsf, NULL); // HDF5_rank = 2
               //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data space created" << std::endl;
               hid_t dataset = H5Dcreate2(unifiedFilePtrSet.h5Var,
                                          "seq_of_vectors",
                                          datatype,
                                          dataspace,
                                          H5P_DEFAULT,  // Link creation property list
                                          H5P_DEFAULT,  // Dataset creation property list
                                          H5P_DEFAULT); // Dataset access property list
               //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data set created" << std::endl;
 
               struct timeval timevalBegin;
               int iRC = UQ_OK_RC;
               iRC = gettimeofday(&timevalBegin,NULL);
               if (iRC) {}; // just to remove compiler warning
 
               //double* dataOut[numParams]; // avoid compiler warning
         std::vector<double*> dataOut((size_t) numParams,NULL);
               dataOut[0] = (double*) malloc(numParams*chainSize*sizeof(double));
               for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1'
                 dataOut[i] = dataOut[i-1] + chainSize; // Yes, just 'chainSize', not 'chainSize*sizeof(double)'
               }
               //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, memory allocated" << std::endl;
               for (unsigned int j = 0; j < chainSize; ++j) {
                 T tmpScalar = m_seq[j];
                 for (unsigned int i = 0; i < numParams; ++i) {
                   dataOut[i][j] = tmpScalar;
                 }
               }
               //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, memory filled" << std::endl;
 
               herr_t status;
               //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 002 \n" << std::endl;
               status = H5Dwrite(dataset,
                                 H5T_NATIVE_DOUBLE,
                                 H5S_ALL,
                                 H5S_ALL,
                                 H5P_DEFAULT,
                                 (void*) dataOut[0]);
               if (status) {}; // just to remove compiler warning
 
               //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 003 \n" << std::endl;
               //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data written" << std::endl;
 
               double writeTime = MiscGetEllapsedSeconds(&timevalBegin);
               if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
                 *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedWriteContents()"
                                         << ": worldRank "      << m_env.worldRank()
                                         << ", fullRank "       << m_env.fullRank()
                                         << ", subEnvironment " << m_env.subId()
                                         << ", subRank "        << m_env.subRank()
                                         << ", inter0Rank "     << m_env.inter0Rank()
                                         << ", fileName = "     << fileName
                                         << ", numParams = "    << numParams
                                         << ", chainSize = "    << chainSize
                                         << ", writeTime = "    << writeTime << " seconds"
                                         << std::endl;
               }
 
               H5Dclose(dataset);
               //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data set closed" << std::endl;
               H5Sclose(dataspace);
               //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data space closed" << std::endl;
               H5Tclose(datatype);
               //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data type closed" << std::endl;
               //free(dataOut[0]); // related to the changes above for compiler warning
               for (unsigned int tmpIndex = 0; tmpIndex < dataOut.size(); tmpIndex++) {
                 free (dataOut[tmpIndex]);
               }
             }
             else {
               queso_error_msg("hdf file type not supported for multiple sub-environments yet");
             }
           }
 #endif
         } // if (m_env.openUnifiedOutputFile())
         //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 004 \n" << std::endl;
       } // if (m_env.inter0Rank() == (int) r)
       m_env.inter0Comm().Barrier();
     } // for r
 
     if (m_env.inter0Rank() == 0) {
       if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
         FilePtrSetStruct unifiedFilePtrSet;
         if (m_env.openUnifiedOutputFile(fileName,
                                         fileType,
                                         false, // Yes, 'writeOver = false' in order to close the array for matlab
                                         unifiedFilePtrSet)) {
           *unifiedFilePtrSet.ofsVar << "];\n";
           m_env.closeFile(unifiedFilePtrSet,fileType);
         }
       }
       else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
         // Do nothing
       }
       else {
         queso_error_msg("invalid file type");
       }
     }
   } // if (m_env.inter0Rank() >= 0)
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
     *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedWriteContents()"
                             << ", fileName = " << fileName
                             << ", fileType = " << fileType
                             << std::endl;
   }
   //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
 
   return;
 }

Member Data Documentation

template<class T = double>

const BaseEnvironment& QUESO::ScalarSequence< T >::m_env

private

Definition at line 516 of file ScalarSequence.h.

template<class T = double>

std::string QUESO::ScalarSequence< T >::m_name

private

Definition at line 517 of file ScalarSequence.h.

Referenced by QUESO::ScalarSequence< T >::copy().

template<class T = double>

std::vector<T> QUESO::ScalarSequence< T >::m_seq

private

Definition at line 518 of file ScalarSequence.h.

Referenced by QUESO::ScalarSequence< T >::append(), and QUESO::ScalarSequence< T >::copy().

template<class T = double>

T* QUESO::ScalarSequence< T >::m_subMaxPlain

mutableprivate

Definition at line 522 of file ScalarSequence.h.

template<class T = double>

T* QUESO::ScalarSequence< T >::m_subMeanPlain

mutableprivate

Definition at line 524 of file ScalarSequence.h.

template<class T = double>

T* QUESO::ScalarSequence< T >::m_subMedianPlain

mutableprivate

Definition at line 526 of file ScalarSequence.h.

template<class T = double>

T* QUESO::ScalarSequence< T >::m_subMinPlain

mutableprivate

Definition at line 520 of file ScalarSequence.h.

template<class T = double>

T* QUESO::ScalarSequence< T >::m_subSampleVariancePlain

mutableprivate

Definition at line 528 of file ScalarSequence.h.

template<class T = double>

T* QUESO::ScalarSequence< T >::m_unifiedMaxPlain

mutableprivate

Definition at line 523 of file ScalarSequence.h.

template<class T = double>

T* QUESO::ScalarSequence< T >::m_unifiedMeanPlain

mutableprivate

Definition at line 525 of file ScalarSequence.h.

template<class T = double>

T* QUESO::ScalarSequence< T >::m_unifiedMedianPlain

mutableprivate

Definition at line 527 of file ScalarSequence.h.

template<class T = double>

T* QUESO::ScalarSequence< T >::m_unifiedMinPlain

mutableprivate

Definition at line 521 of file ScalarSequence.h.

template<class T = double>

T* QUESO::ScalarSequence< T >::m_unifiedSampleVariancePlain

mutableprivate

Definition at line 529 of file ScalarSequence.h.

The documentation for this class was generated from the following files:

src/basic/inc/ScalarSequence.h
src/basic/src/ScalarSequence.C

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

template<class T = double> class QUESO::ScalarSequence< T >

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

template<class T = double>
class QUESO::ScalarSequence< T >