queso-0.56.0
Private Member Functions | Private Attributes | List of all members
QUESO::ScalarSequence< T > Class Template Reference

Class for handling scalar samples. More...

#include <ScalarSequence.h>

Collaboration diagram for QUESO::ScalarSequence< T >:
Collaboration graph
[legend]

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 BaseEnvironmentenv () 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...
 
subMeanExtra (unsigned int initialPos, unsigned int numPos) const
 Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPos. More...
 
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...
 
subMedianExtra (unsigned int initialPos, unsigned int numPos) const
 Finds the median value of the sub-sequence, considering numPos positions starting at position initialPos. More...
 
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...
 
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...
 
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...
 
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...
 
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...
 
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...
 
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...
 
autoCovariance (unsigned int initialPos, unsigned int numPos, const T &meanValue, unsigned int lag) const
 Calculates the autocovariance. More...
 
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...
 
subInterQuantileRange (unsigned int initialPos) const
 Returns the interquartile range of the values in the sub-sequence. More...
 
unifiedInterQuantileRange (bool useOnlyInter0Comm, unsigned int initialPos) const
 Returns the interquartile range of the values in the unified sequence. More...
 
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...
 
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...
 
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...
 
subPositionsOfMaximum (const ScalarSequence< T > &subCorrespondingScalarValues, ScalarSequence< T > &subPositionsOfMaximum)
 Finds the positions where the maximum element occurs in the sub-sequence. More...
 
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 writeUnifiedMatlabHeader (std::ofstream &ofs, double sequenceSize) const
 Helper function to write header info for matlab files from all chains. More...
 
void writeSubMatlabHeader (std::ofstream &ofs, double sequenceSize) const
 Helper function to write header info for matlab files from one chain. More...
 
void writeTxtHeader (std::ofstream &ofs, double sequenceSize) const
 Helper function to write txt info for matlab files. 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 BaseEnvironmentm_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.

36  :
37  m_env (env),
38  m_name (name),
40  m_subMinPlain (NULL),
41  m_unifiedMinPlain (NULL),
42  m_subMaxPlain (NULL),
43  m_unifiedMaxPlain (NULL),
44  m_subMeanPlain (NULL),
45  m_unifiedMeanPlain (NULL),
46  m_subMedianPlain (NULL),
47  m_unifiedMedianPlain (NULL),
50 {
51 }
const std::string & name() const
Access to the name of the sequence of scalars.
std::vector< T > m_seq
const BaseEnvironment & env() const
Access to QUESO environment.
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
template<class T >
QUESO::ScalarSequence< T >::~ScalarSequence ( )

Destructor.

Definition at line 54 of file ScalarSequence.C.

55 {
57 }
void deleteStoredScalars()
Deletes all stored scalars.

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 2488 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().

2492 {
2493  queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+1), "srcInitialPos is too big");
2494 
2495  queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+srcNumPos), "srcNumPos is too big");
2496 
2498  unsigned int currentSize = this->subSequenceSize();
2499  m_seq.resize(currentSize+srcNumPos,0.);
2500  for (unsigned int i = 0; i < srcNumPos; ++i) {
2501  m_seq[currentSize+i] = src.m_seq[srcInitialPos+i];
2502  }
2503 
2504  return;
2505 }
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
void deleteStoredScalars()
Deletes all stored scalars.
std::vector< T > m_seq
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1334 of file ScalarSequence.C.

References queso_require_msg.

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

1338 {
1339  bool bRC = ((initialPos < this->subSequenceSize()) &&
1340  (0 < numPos ) &&
1341  ((initialPos+numPos) <= this->subSequenceSize()) &&
1342  (lag < numPos )); // lag should not be too large
1343  queso_require_msg(bRC, "invalid input data");
1344 
1345  T meanValue = this->subMeanExtra(initialPos,
1346  numPos);
1347 
1348  T covValueZero = this->autoCovariance(initialPos,
1349  numPos,
1350  meanValue,
1351  0); // lag
1352 
1353  T corrValue = this->autoCovariance(initialPos,
1354  numPos,
1355  meanValue,
1356  lag);
1357 
1358  return corrValue/covValueZero;
1359 }
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
T autoCovariance(unsigned int initialPos, unsigned int numPos, const T &meanValue, unsigned int lag) const
Calculates the autocovariance.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1363 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().

1368 {
1369  double tmp = log((double) numPos)/log(2.);
1370  double fractionalPart = tmp - ((double) ((unsigned int) tmp));
1371  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
1372  unsigned int fftSize = (unsigned int) std::pow(2.,tmp+1);
1373 
1374  std::vector<double> rawDataVec(numPos,0.);
1375  std::vector<std::complex<double> > resultData(0,std::complex<double>(0.,0.));
1376  Fft<T> fftObj(m_env);
1377 
1378  // Forward FFT
1379  this->extractRawData(initialPos,
1380  1, // spacing
1381  numPos,
1382  rawDataVec);
1383  T meanValue = this->subMeanExtra(initialPos,
1384  numPos);
1385  for (unsigned int j = 0; j < numPos; ++j) {
1386  rawDataVec[j] -= meanValue; // IMPORTANT
1387  }
1388 
1389  rawDataVec.resize(fftSize,0.);
1390 
1391  //if (m_env.subDisplayFile()) {
1392  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1393  // << ": about to call fftObj.forward()"
1394  // << " with rawDataVec.size() = " << rawDataVec.size()
1395  // << ", fftSize = " << fftSize
1396  // << ", resultData.size() = " << resultData.size()
1397  // << std::endl;
1398  //}
1399  fftObj.forward(rawDataVec,fftSize,resultData);
1400 
1401  // Inverse FFT
1402  for (unsigned int j = 0; j < fftSize; ++j) {
1403  rawDataVec[j] = std::norm(resultData[j]);
1404  }
1405  //if (m_env.subDisplayFile()) {
1406  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1407  // << ": about to call fftObj.inverse()"
1408  // << " with rawDataVec.size() = " << rawDataVec.size()
1409  // << ", fftSize = " << fftSize
1410  // << ", resultData.size() = " << resultData.size()
1411  // << std::endl;
1412  //}
1413  fftObj.inverse(rawDataVec,fftSize,resultData);
1414  //if (m_env.subDisplayFile()) {
1415  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1416  // << ": returned succesfully from fftObj.inverse()"
1417  // << std::endl;
1418  //}
1419 
1420  // Prepare return data
1421  autoCorrs.resize(maxLag+1,0.); // Yes, +1
1422  for (unsigned int j = 0; j < autoCorrs.size(); ++j) {
1423  double ratio = ((double) j)/((double) (numPos-1));
1424  autoCorrs[j] = ( resultData[j].real()/resultData[0].real() )*(1.-ratio);
1425  }
1426 
1427  return;
1428 }
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, std::vector< double > &rawData) const
Extracts the raw data.
const BaseEnvironment & m_env
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 1432 of file ScalarSequence.C.

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

1437 {
1438  //if (m_env.subDisplayFile()) {
1439  // *m_env.subDisplayFile() << "Entering ScalarSequence<T>::autoCorrViaFft(), for sum"
1440  // << ": initialPos = " << initialPos
1441  // << ", numPos = " << numPos
1442  // << std::endl;
1443  //}
1444 
1445  double tmp = log((double) numPos)/log(2.);
1446  double fractionalPart = tmp - ((double) ((unsigned int) tmp));
1447  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
1448  unsigned int fftSize = (unsigned int) std::pow(2.,tmp+1);
1449 
1450  std::vector<double> rawDataVec(numPos,0.);
1451  std::vector<std::complex<double> > resultData(0,std::complex<double>(0.,0.));
1452  Fft<T> fftObj(m_env);
1453 
1454  // Forward FFT
1455  this->extractRawData(initialPos,
1456  1, // spacing
1457  numPos,
1458  rawDataVec);
1459  T meanValue = this->subMeanExtra(initialPos,
1460  numPos);
1461  for (unsigned int j = 0; j < numPos; ++j) {
1462  rawDataVec[j] -= meanValue; // IMPORTANT
1463  }
1464  rawDataVec.resize(fftSize,0.);
1465 
1466  //if (m_env.subDisplayFile()) {
1467  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft(), for sum"
1468  // << ": about to call fftObj.forward()"
1469  // << " with rawDataVec.size() = " << rawDataVec.size()
1470  // << ", fftSize = " << fftSize
1471  // << ", resultData.size() = " << resultData.size()
1472  // << std::endl;
1473  //}
1474  fftObj.forward(rawDataVec,fftSize,resultData);
1475 
1476  // Inverse FFT
1477  for (unsigned int j = 0; j < fftSize; ++j) {
1478  rawDataVec[j] = std::norm(resultData[j]);
1479  }
1480  fftObj.inverse(rawDataVec,fftSize,resultData);
1481 
1482  //if (m_env.subDisplayFile()) {
1483  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft(), for sum"
1484  // << ": computed auto covariance for lag 0 = " << resultData[0].real()/((double) (numPos))
1485  // << ", computed resultData[0].imag() = " << resultData[0].imag()
1486  // << std::endl;
1487  //}
1488 
1489  // Prepare return data
1490  autoCorrsSum = 0.;
1491  for (unsigned int j = 0; j < numSum; ++j) { // Yes, begin at lag '0'
1492  double ratio = ((double) j)/((double) (numPos-1));
1493  autoCorrsSum += ( resultData[j].real()/resultData[0].real() )*(1.-ratio);
1494  }
1495 
1496  return;
1497 }
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, std::vector< double > &rawData) const
Extracts the raw data.
const BaseEnvironment & m_env
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 1304 of file ScalarSequence.C.

References queso_require_msg.

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

1309 {
1310  bool bRC = ((initialPos < this->subSequenceSize()) &&
1311  (0 < numPos ) &&
1312  ((initialPos+numPos) <= this->subSequenceSize()) &&
1313  (lag < numPos )); // lag should not be too large
1314  queso_require_msg(bRC, "invalid input data");
1315 
1316  unsigned int loopSize = numPos - lag;
1317  unsigned int finalPosPlus1 = initialPos + loopSize;
1318  T diff1;
1319  T diff2;
1320  T covValue = 0.;
1321  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1322  diff1 = m_seq[j ] - meanValue;
1323  diff2 = m_seq[j+lag] - meanValue;
1324  covValue += diff1*diff2;
1325  }
1326 
1327  covValue /= (T) loopSize;
1328 
1329  return covValue;
1330 }
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2459 of file ScalarSequence.C.

References queso_error_msg, and queso_not_implemented.

2463 {
2464  double resultValue = 0.;
2465 
2466  // FIX ME: put logic if (numSubEnvs == 1) ...
2467 
2468  if (useOnlyInter0Comm) {
2469  if (m_env.inter0Rank() >= 0) {
2471  }
2472  else {
2473  // Node not in the 'inter0' communicator
2474  // Do nothing
2475  }
2476  }
2477  else {
2478  queso_error_msg("parallel vectors not supported yet");
2479  }
2480 
2481  //m_env.fullComm().Barrier();
2482 
2483  return resultValue;
2484 }
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
#define queso_not_implemented()
Definition: asserts.h:56
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
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().

124 {
125  unsigned int numPos = this->subSequenceSize();
126  if (numPos) {
127  this->resetValues(0,numPos);
128  this->resizeSequence(0);
129  }
130 
131  return;
132 }
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
void resetValues(unsigned int initialPos, unsigned int)
Sets numPos values of the sequence to zero, starting at position initialPos.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 3168 of file ScalarSequence.C.

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

3169 {
3170  m_name = src.m_name;
3171  m_seq.clear();
3172  m_seq.resize(src.subSequenceSize(),0.);
3173  for (unsigned int i = 0; i < m_seq.size(); ++i) {
3174  m_seq[i] = src.m_seq[i];
3175  }
3177 
3178  return;
3179 }
void deleteStoredScalars()
Deletes all stored scalars.
std::vector< T > m_seq
template<class T >
void QUESO::ScalarSequence< T >::deleteStoredScalars ( )

Deletes all stored scalars.

Definition at line 467 of file ScalarSequence.C.

468 {
469  if (m_subMinPlain) {
470  delete m_subMinPlain;
471  m_subMinPlain = NULL;
472  }
473  if (m_unifiedMinPlain) {
474  delete m_unifiedMinPlain;
475  m_unifiedMinPlain = NULL;
476  }
477  if (m_subMaxPlain) {
478  delete m_subMaxPlain;
479  m_subMaxPlain = NULL;
480  }
481  if (m_unifiedMaxPlain) {
482  delete m_unifiedMaxPlain;
483  m_unifiedMaxPlain = NULL;
484  }
485  if (m_subMeanPlain) {
486  delete m_subMeanPlain;
487  m_subMeanPlain = NULL;
488  }
489  if (m_unifiedMeanPlain) {
490  delete m_unifiedMeanPlain;
491  m_unifiedMeanPlain = NULL;
492  }
493  if (m_subMedianPlain) {
494  delete m_subMedianPlain;
495  m_subMedianPlain = NULL;
496  }
497  if (m_unifiedMedianPlain) {
498  delete m_unifiedMedianPlain;
499  m_unifiedMedianPlain = NULL;
500  }
504  }
508  }
509 
510  return;
511 }
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().

102 {
103  return m_env;
104 }
const BaseEnvironment & 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().

207 {
208  if (this->subSequenceSize() == 0) return;
209 
210  bool bRC = ((initialPos < this->subSequenceSize()) &&
211  (0 < numPos ) &&
212  ((initialPos+numPos) <= this->subSequenceSize()));
213  queso_require_msg(bRC, "invalid input data");
214 
215  seqScalarPositionIteratorTypedef posIteratorBegin = m_seq.begin();
216  if (initialPos < this->subSequenceSize()) std::advance(posIteratorBegin,initialPos);
217  else posIteratorBegin = m_seq.end();
218 
219  unsigned int posEnd = initialPos + numPos - 1;
220  seqScalarPositionIteratorTypedef posIteratorEnd = m_seq.begin();
221  if (posEnd < this->subSequenceSize()) std::advance(posIteratorEnd,posEnd);
222  else posIteratorEnd = m_seq.end();
223 
224  unsigned int oldSequenceSize = this->subSequenceSize();
225  m_seq.erase(posIteratorBegin,posIteratorEnd);
226  queso_require_equal_to_msg((oldSequenceSize - numPos), this->subSequenceSize(), "(oldSequenceSize - numPos) != this->subSequenceSize()");
227 
229 
230  return;
231 }
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
void deleteStoredScalars()
Deletes all stored scalars.
std::vector< T > m_seq
std::vector< T >::iterator seqScalarPositionIteratorTypedef
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 3224 of file ScalarSequence.C.

3229 {
3230  rawDataVec.resize(numPos);
3231  if (spacing == 1) {
3232  for (unsigned int j = 0; j < numPos; ++j) {
3233  rawDataVec[j] = m_seq[initialPos+j ];
3234  }
3235  }
3236  else {
3237  for (unsigned int j = 0; j < numPos; ++j) {
3238  rawDataVec[j] = m_seq[initialPos+j*spacing];
3239  }
3240  }
3241 
3242  return;
3243 }
std::vector< T > m_seq
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 3183 of file ScalarSequence.C.

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

3188 {
3189  scalarSeq.resizeSequence(numPos);
3190  if (spacing == 1) {
3191  for (unsigned int j = 0; j < numPos; ++j) {
3192  //if ((initialPos+j*spacing) > m_seq.size()) {
3193  // std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
3194  // << ": initialPos = " << initialPos
3195  // << ", spacing = " << spacing
3196  // << ", numPos = " << numPos
3197  // << ", j = " << j
3198  // << ", position got too large"
3199  // << std::endl;
3200  //}
3201  scalarSeq[j] = m_seq[initialPos+j ];
3202  }
3203  }
3204  else {
3205  for (unsigned int j = 0; j < numPos; ++j) {
3206  //if ((initialPos+j*spacing) > m_seq.size()) {
3207  // std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
3208  // << ": initialPos = " << initialPos
3209  // << ", spacing = " << spacing
3210  // << ", numPos = " << numPos
3211  // << ", j = " << j
3212  // << ", position got too large"
3213  // << std::endl;
3214  //}
3215  scalarSeq[j] = m_seq[initialPos+j*spacing];
3216  }
3217  }
3218 
3219  return;
3220 }
std::vector< T > m_seq
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 2420 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().

2423 {
2424  if (m_env.subDisplayFile()) {
2425  *m_env.subDisplayFile() << "Entering ScalarSequence<V,M>::filter()"
2426  << ": initialPos = " << initialPos
2427  << ", spacing = " << spacing
2428  << ", subSequenceSize = " << this->subSequenceSize()
2429  << std::endl;
2430  }
2431 
2432  unsigned int i = 0;
2433  unsigned int j = initialPos;
2434  unsigned int originalSubSequenceSize = this->subSequenceSize();
2435  while (j < originalSubSequenceSize) {
2436  if (i != j) {
2437  //*m_env.subDisplayFile() << i << "--" << j << " ";
2438  m_seq[i] = m_seq[j];
2439  }
2440  i++;
2441  j += spacing;
2442  }
2443 
2444  this->resizeSequence(i);
2445 
2446  if (m_env.subDisplayFile()) {
2447  *m_env.subDisplayFile() << "Leaving ScalarSequence<V,M>::filter()"
2448  << ": initialPos = " << initialPos
2449  << ", spacing = " << spacing
2450  << ", subSequenceSize = " << this->subSequenceSize()
2451  << std::endl;
2452  }
2453 
2454  return;
2455 }
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
std::vector< T > m_seq
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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, and RawValue_MPI_DOUBLE.

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

238 {
239  // The logic (numSubEnvs == 1) does *not* apply here because 'outputVec' needs to be filled
240  //if (m_env.numSubEnvironments() == 1) {
241  // // No need to do anything
242  // return;
243  //}
244 
245  if (useOnlyInter0Comm) {
246  if (m_env.inter0Rank() >= 0) {
247  int auxSubSize = (int) this->subSequenceSize();
248  unsigned int auxUnifiedSize = this->unifiedSequenceSize(useOnlyInter0Comm);
249  outputVec.resize(auxUnifiedSize,0.);
250 
251  //******************************************************************
252  // Use MPI_Gatherv for the case different nodes have different amount of data // KAUST4
253  //******************************************************************
254  std::vector<int> recvcnts(m_env.inter0Comm().NumProc(),0); // '0' is NOT the correct value for recvcnts[0]
255  m_env.inter0Comm().template Gather<int>(&auxSubSize, 1, &recvcnts[0], (int) 1, 0,
256  "ScalarSequence<T>::getUnifiedContentsAtProc0Only()",
257  "failed MPI.Gather()");
258  if (m_env.inter0Rank() == 0) {
259  //recvcnts[0] = (int) this->subSequenceSize(); // FIX ME: really necessary????
260  queso_require_equal_to_msg(recvcnts[0], (int) this->subSequenceSize(), "failed MPI.Gather() result at proc 0");
261  }
262 
263  std::vector<int> displs(m_env.inter0Comm().NumProc(),0);
264  for (unsigned int r = 1; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) { // Yes, from '1' on
265  displs[r] = displs[r-1] + recvcnts[r-1];
266  }
267 
268 #if 0 // for debug only
269  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
270  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) {
271  *m_env.subDisplayFile() << " auxSubSize = " << auxSubSize
272  << ", recvcnts[" << r << "] = " << recvcnts[r]
273  << ", displs[" << r << "] = " << displs[r]
274  << ", m_seq.size() = " << m_seq.size()
275  << ", outputVec.size() = " << outputVec.size()
276  << std::endl;
277  }
278  for (unsigned int i = 0; i < m_seq.size(); ++i) {
279  *m_env.subDisplayFile() << " (before gatherv) m_seq[" << i << "]= " << m_seq[i]
280  << std::endl;
281  }
282  }
283 #endif
284 
285  m_env.inter0Comm().template Gatherv<double>(&m_seq[0], auxSubSize,
286  &outputVec[0], (int *) &recvcnts[0], (int *) &displs[0], 0,
287  "ScalarSequence<T>::getUnifiedContentsAtProc0Only()",
288  "failed MPI.Gatherv()");
289 
290 #if 0 // for debug only
291  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
292  for (unsigned int i = 0; i < m_seq.size(); ++i) {
293  *m_env.subDisplayFile() << " (after gatherv) m_seq[" << i << "]= " << m_seq[i]
294  << std::endl;
295  }
296  for (unsigned int i = 0; i < outputVec.size(); ++i) {
297  *m_env.subDisplayFile() << " (after gatherv) outputVec[" << i << "]= " << outputVec[i]
298  << std::endl;
299  }
300  }
301 #endif
302 
303 #if 0 // for debug only
304  if (m_env.inter0Rank() == 0) {
305  for (unsigned int i = 0; i < auxSubSize; ++i) {
306  outputVec[i] = m_seq[i];
307  }
308  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,
309  "ScalarSequence<T>::getUnifiedContentsAtProc0Only(1)",
310  "failed MPI.Gatherv()");
311  }
312  else {
313  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,
314  "ScalarSequence<T>::getUnifiedContentsAtProc0Only(2)",
315  "failed MPI.Gatherv()");
316  }
317  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
318  for (unsigned int i = 0; i < m_seq.size(); ++i) {
319  *m_env.subDisplayFile() << " (after 2nd gatherv) m_seq[" << i << "]= " << m_seq[i]
320  << std::endl;
321  }
322  for (unsigned int i = 0; i < outputVec.size(); ++i) {
323  *m_env.subDisplayFile() << " (after 2nd gatherv) outputVec[" << i << "]= " << outputVec[i]
324  << std::endl;
325  }
326  }
327 #endif
328  }
329  else {
330  // Node not in the 'inter0' communicator
331  // No need to do anything
332  }
333  }
334  else {
335  queso_error_msg("parallel vectors not supported yet");
336  }
337 
338  return;
339 }
int NumProc() const
Returns total number of processes.
Definition: MpiComm.C:123
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
void Gatherv(void *sendbuf, int sendcnt, RawType_MPI_Datatype sendtype, void *recvbuf, int *recvcnts, int *displs, RawType_MPI_Datatype recvtype, int root, const char *whereMsg, const char *whatMsg) const
Gathers into specified locations from all processes in a group.
Definition: MpiComm.C:244
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:67
unsigned int displayVerbosity() const
Definition: Environment.C:449
std::vector< T > m_seq
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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.

109 {
110  return m_name;
111 }
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.

62 {
63  this->copy(rhs);
64  return *this;
65 }
void copy(const ScalarSequence< T > &src)
Copies the scalar sequence src to 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.

70 {
71  if (posId >= this->subSequenceSize()) {
72  std::cerr << "In ScalarSequence<T>::operator[]() const"
73  << ": posId = " << posId
74  << ", this->subSequenceSize() = " << this->subSequenceSize()
75  << std::endl;
76  }
77  queso_require_less_msg(posId, this->subSequenceSize(), "posId > subSequenceSize()");
78 
79  return m_seq[posId];
80 }
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:75
std::vector< T > m_seq
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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.

85 {
86  if (posId >= this->subSequenceSize()) {
87  std::cerr << "In ScalarSequence<T>::operator[]()"
88  << ": posId = " << posId
89  << ", this->subSequenceSize() = " << this->subSequenceSize()
90  << std::endl;
91  }
92  queso_require_less_msg(posId, this->subSequenceSize(), "posId > subSequenceSize()");
93 
95 
96  return m_seq[posId];
97 }
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:75
void deleteStoredScalars()
Deletes all stored scalars.
std::vector< T > m_seq
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 3265 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.

3269 {
3270  int parentNode = m_env.inter0Rank() & ~(1 << currentTreeLevel);
3271 
3272  if (m_env.inter0Rank() >= 0) { // KAUST
3273  if (currentTreeLevel == 0) {
3274  // Leaf node: sort own local data.
3275  unsigned int leafDataSize = leafData.size();
3276  sortedBuffer.resize(leafDataSize,0.);
3277  for (unsigned int i = 0; i < leafDataSize; ++i) {
3278  sortedBuffer[i] = leafData[i];
3279  }
3280  std::sort(sortedBuffer.begin(), sortedBuffer.end());
3281  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3282  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3283  << ": tree node " << m_env.inter0Rank()
3284  << ", leaf sortedBuffer[0] = " << sortedBuffer[0]
3285  << ", leaf sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
3286  << std::endl;
3287  }
3288  }
3289  else {
3290  int nextTreeLevel = currentTreeLevel - 1;
3291  int rightChildNode = m_env.inter0Rank() | (1 << nextTreeLevel);
3292 
3293  if (rightChildNode >= m_env.inter0Comm().NumProc()) { // No right child. Move down one level.
3294  this->parallelMerge(sortedBuffer,
3295  leafData,
3296  nextTreeLevel);
3297  }
3298  else {
3299  unsigned int uintBuffer[1];
3300  uintBuffer[0] = nextTreeLevel;
3301  m_env.inter0Comm().Send((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, rightChildNode, SCALAR_SEQUENCE_INIT_MPI_MSG,
3302  "ScalarSequence<T>::parallelMerge()",
3303  "failed MPI.Send() for init");
3304 
3305  this->parallelMerge(sortedBuffer,
3306  leafData,
3307  nextTreeLevel);
3308 
3309  // Prepare variable 'leftSortedBuffer': just copy own current sorted data.
3310  unsigned int leftSize = sortedBuffer.size();
3311  std::vector<T> leftSortedBuffer(leftSize,0.);
3312  for (unsigned int i = 0; i < leftSize; ++i) {
3313  leftSortedBuffer[i] = sortedBuffer[i];
3314  }
3315 
3316  // Prepare variable 'rightSortedBuffer': receive data from right child node.
3317  RawType_MPI_Status status;
3318  m_env.inter0Comm().Recv((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, rightChildNode, SCALAR_SEQUENCE_SIZE_MPI_MSG, &status,
3319  "ScalarSequence<T>::parallelMerge()",
3320  "failed MPI.Recv() for size");
3321  //if (status) {}; // just to remove compiler warning
3322 
3323  unsigned int rightSize = uintBuffer[0];
3324  std::vector<T> rightSortedBuffer(rightSize,0.);
3325  m_env.inter0Comm().Recv((void *) &rightSortedBuffer[0], (int) rightSize, RawValue_MPI_DOUBLE, rightChildNode, SCALAR_SEQUENCE_DATA_MPI_MSG, &status,
3326  "ScalarSequence<T>::parallelMerge()",
3327  "failed MPI.Recv() for data");
3328 
3329  // Merge the two results into 'sortedBuffer'.
3330  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3331  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3332  << ": tree node " << m_env.inter0Rank()
3333  << " is combining " << leftSortedBuffer.size()
3334  << " left doubles with " << rightSortedBuffer.size()
3335  << " right doubles"
3336  << std::endl;
3337  }
3338 
3339  sortedBuffer.clear();
3340  sortedBuffer.resize(leftSortedBuffer.size()+rightSortedBuffer.size(),0.);
3341  unsigned int i = 0;
3342  unsigned int j = 0;
3343  unsigned int k = 0;
3344  while ((i < leftSize ) &&
3345  (j < rightSize)) {
3346  if (leftSortedBuffer[i] > rightSortedBuffer[j]) sortedBuffer[k++] = rightSortedBuffer[j++];
3347  else sortedBuffer[k++] = leftSortedBuffer [i++];
3348  }
3349  while (i < leftSize ) sortedBuffer[k++] = leftSortedBuffer [i++];
3350  while (j < rightSize) sortedBuffer[k++] = rightSortedBuffer[j++];
3351  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3352  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3353  << ": tree node " << m_env.inter0Rank()
3354  << ", merged sortedBuffer[0] = " << sortedBuffer[0]
3355  << ", merged sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
3356  << std::endl;
3357  }
3358  }
3359  }
3360 
3361  if (parentNode != m_env.inter0Rank()) {
3362  // Transmit data to parent node.
3363  unsigned int uintBuffer[1];
3364  uintBuffer[0] = sortedBuffer.size();
3365  m_env.inter0Comm().Send((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, parentNode, SCALAR_SEQUENCE_SIZE_MPI_MSG,
3366  "ScalarSequence<T>::parallelMerge()",
3367  "failed MPI.Send() for size");
3368 
3369  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
3370  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3371  << ": tree node " << m_env.inter0Rank()
3372  << " is sending " << sortedBuffer.size()
3373  << " doubles to tree node " << parentNode
3374  << ", with sortedBuffer[0] = " << sortedBuffer[0]
3375  << " and sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
3376  << std::endl;
3377  }
3378 
3379  m_env.inter0Comm().Send((void *) &sortedBuffer[0], (int) sortedBuffer.size(), RawValue_MPI_DOUBLE, parentNode, SCALAR_SEQUENCE_DATA_MPI_MSG,
3380  "ScalarSequence<T>::parallelMerge()",
3381  "failed MPI.Send() for data");
3382  }
3383  } // KAUST
3384 
3385  return;
3386 }
int NumProc() const
Returns total number of processes.
Definition: MpiComm.C:123
int RawType_MPI_Status
Definition: MpiComm.h:62
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
void Send(void *buf, int count, RawType_MPI_Datatype datatype, int dest, int tag, const char *whereMsg, const char *whatMsg) const
Possibly blocking send of data from this process to another process.
Definition: MpiComm.C:310
void parallelMerge(std::vector< T > &sortedBuffer, const std::vector< T > &leafData, unsigned int treeLevel) const
Sorts/merges data in parallel using MPI.
int k
Definition: ann_sample.cpp:53
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:67
#define RawValue_MPI_UNSIGNED
Definition: MpiComm.h:68
#define SCALAR_SEQUENCE_INIT_MPI_MSG
unsigned int displayVerbosity() const
Definition: Environment.C:449
void Recv(void *buf, int count, RawType_MPI_Datatype datatype, int source, int tag, RawType_MPI_Status *status, const char *whereMsg, const char *whatMsg) const
Blocking receive of data from this process to another process.
Definition: MpiComm.C:297
#define SCALAR_SEQUENCE_SIZE_MPI_MSG
const BaseEnvironment & m_env
#define SCALAR_SEQUENCE_DATA_MPI_MSG
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
template<class T >
std::vector< T > & QUESO::ScalarSequence< T >::rawData ( )
private

The sequence of scalars. Access to private attribute m_seq.

Definition at line 3247 of file ScalarSequence.C.

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

3248 {
3249  return m_seq;
3250 }
std::vector< T > 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.

187 {
188  if (this->subSequenceSize() == 0) return;
189 
190  bool bRC = ((initialPos < this->subSequenceSize()) &&
191  (0 < numPos ) &&
192  ((initialPos+numPos) <= this->subSequenceSize()));
193  queso_require_msg(bRC, "invalid input data");
194 
195  for (unsigned int j = 0; j < numPos; ++j) {
196  m_seq[initialPos+j] = 0.;
197  }
198 
200 
201  return;
202 }
void deleteStoredScalars()
Deletes all stored scalars.
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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().

174 {
175  if (newSequenceSize != this->subSequenceSize()) {
176  m_seq.resize(newSequenceSize,0.);
177  std::vector<T>(m_seq).swap(m_seq);
179  }
180 
181  return;
182 }
void deleteStoredScalars()
Deletes all stored scalars.
std::vector< T > m_seq
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 515 of file ScalarSequence.C.

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

516 {
517  unsigned int maxJ = this->subSequenceSize();
518  if (meanValue == 0.) {
519  for (unsigned int j = 0; j < maxJ; ++j) {
520  m_seq[j] = m_env.rngObject()->gaussianSample(stdDev);
521  }
522  }
523  else {
524  for (unsigned int j = 0; j < maxJ; ++j) {
525  m_seq[j] = meanValue + m_env.rngObject()->gaussianSample(stdDev);
526  }
527  }
528 
530 
531  return;
532 }
const RngBase * rngObject() const
Access to the RNG object.
Definition: Environment.C:470
void deleteStoredScalars()
Deletes all stored scalars.
std::vector< T > m_seq
virtual double gaussianSample(double stdDev) const =0
Samples a value from a Gaussian distribution with standard deviation given by stdDev.
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
template<class T >
void QUESO::ScalarSequence< T >::setName ( const std::string &  newName)
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 536 of file ScalarSequence.C.

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

537 {
538  unsigned int maxJ = this->subSequenceSize();
539  if (a == 0.) {
540  if (b == 1.) {
541  for (unsigned int j = 0; j < maxJ; ++j) {
543  }
544  }
545  else {
546  for (unsigned int j = 0; j < maxJ; ++j) {
547  m_seq[j] = b*m_env.rngObject()->uniformSample();
548  }
549  }
550  }
551  else {
552  if ((b-a) == 1.) {
553  for (unsigned int j = 0; j < maxJ; ++j) {
554  m_seq[j] = a + m_env.rngObject()->uniformSample();
555  }
556  }
557  else {
558  for (unsigned int j = 0; j < maxJ; ++j) {
559  m_seq[j] = a + (b-a)*m_env.rngObject()->uniformSample();
560  }
561  }
562  }
563 
565 
566  return;
567 }
const RngBase * rngObject() const
Access to the RNG object.
Definition: Environment.C:470
virtual double uniformSample() const =0
Samples a value from a uniform distribution.
void deleteStoredScalars()
Deletes all stored scalars.
std::vector< T > m_seq
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 712 of file ScalarSequence.C.

716 {
717  T tmpMinValue;
718  T tmpMaxValue;
719  std::vector<unsigned int> bins(numEvaluationPoints,0);
720 
721  subMinMaxExtra(0, // initialPos
722  this->subSequenceSize(),
723  tmpMinValue,
724  tmpMaxValue);
725  subBasicHistogram(0, // initialPos,
726  tmpMinValue,
727  tmpMaxValue,
728  gridValues,
729  bins);
730 
731  unsigned int sumOfBins = 0;
732  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
733  sumOfBins += bins[i];
734  }
735 
736  cdfValues.clear();
737  cdfValues.resize(numEvaluationPoints);
738  unsigned int partialSum = 0;
739  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
740  partialSum += bins[i];
741  cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
742  }
743 
744  return;
745 }
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.
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1738 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

1744 {
1745  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1746 
1747  for (unsigned int j = 0; j < bins.size(); ++j) {
1748  bins[j] = 0;
1749  }
1750 
1751  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1752  double minCenter = minHorizontalValue - horizontalDelta/2.;
1753  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1754  gridValues = new UniformOneDGrid<T>(m_env,
1755  "",
1756  bins.size(),
1757  minCenter,
1758  maxCenter);
1759 
1760  unsigned int dataSize = this->subSequenceSize();
1761  for (unsigned int j = 0; j < dataSize; ++j) {
1762  double value = m_seq[j];
1763  if (value < minHorizontalValue) {
1764  bins[0] += value;
1765  }
1766  else if (value >= maxHorizontalValue) {
1767  bins[bins.size()-1] += value;
1768  }
1769  else {
1770  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1771  bins[index] += value;
1772  }
1773  }
1774 
1775  return;
1776 }
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
std::vector< T > m_seq
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2310 of file ScalarSequence.C.

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

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

2315 {
2316  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2317  (0 < evaluationPositions.size()) &&
2318  (evaluationPositions.size() == densityValues.size() ));
2319  queso_require_msg(bRC, "invalid input data");
2320 
2321  unsigned int dataSize = this->subSequenceSize() - initialPos;
2322  unsigned int numEvals = evaluationPositions.size();
2323 
2324  double scaleInv = 1./scaleValue;
2325  for (unsigned int j = 0; j < numEvals; ++j) {
2326  double x = evaluationPositions[j];
2327  double value = 0.;
2328  for (unsigned int k = 0; k < dataSize; ++k) {
2329  double xk = m_seq[initialPos+k];
2330  value += MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
2331  }
2332  densityValues[j] = scaleInv * (value/(double) dataSize);
2333  }
2334 
2335  return;
2336 }
double MiscGaussianDensity(double x, double mu, double sigma)
int k
Definition: ann_sample.cpp:53
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1601 of file ScalarSequence.C.

References queso_require_equal_to_msg, and queso_require_greater_equal_msg.

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

1607 {
1608  queso_require_equal_to_msg(centers.size(), bins.size(), "vectors 'centers' and 'bins' have different sizes");
1609 
1610  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1611 
1612  if (initialPos) {}; // just to remove compiler warning
1613 
1614  for (unsigned int j = 0; j < bins.size(); ++j) {
1615  centers[j] = 0.;
1616  bins[j] = 0;
1617  }
1618 
1619  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1620 
1621  double minCenter = minHorizontalValue - horizontalDelta/2.;
1622  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1623  for (unsigned int j = 0; j < centers.size(); ++j) {
1624  double factor = ((double) j)/(((double) centers.size()) - 1.);
1625  centers[j] = (1. - factor) * minCenter + factor * maxCenter;
1626  }
1627 
1628  unsigned int dataSize = this->subSequenceSize();
1629  for (unsigned int j = 0; j < dataSize; ++j) {
1630  double value = m_seq[j];
1631  if (value < minHorizontalValue) {
1632  bins[0]++;
1633  }
1634  else if (value >= maxHorizontalValue) {
1635  bins[bins.size()-1]++;
1636  }
1637  else {
1638  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1639  bins[index]++;
1640  }
1641  }
1642 
1643  return;
1644 }
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
std::vector< T > m_seq
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1995 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().

1996 {
1997  queso_require_less_msg(initialPos, this->subSequenceSize(), "'initialPos' is too big");
1998 
1999  ScalarSequence sortedSequence(m_env,0,"");
2000  this->subSort(initialPos,
2001  sortedSequence);
2002 
2003  // The test above guarantees that 'dataSize >= 1'
2004  unsigned int dataSize = this->subSequenceSize() - initialPos;
2005 
2006  queso_require_equal_to_msg(dataSize, sortedSequence.subSequenceSize(), "inconsistent size variables");
2007 
2008  bool everythingOk = true;
2009 
2010  // pos1 = (dataSize+1)/4 - 1
2011  // pos1 >= 0 <==> dataSize >= 3
2012  // pos1 < (dataSize-1) <==> 3*dataSize > 1
2013  unsigned int pos1 = (unsigned int) ( (((double) dataSize) + 1.)*1./4. - 1. );
2014  if (pos1 > (dataSize-1)) {
2015  pos1 = 0;
2016  everythingOk = false;
2017  }
2018  unsigned int pos1inc = pos1+1;
2019  if (pos1inc > (dataSize-1)) {
2020  pos1inc = dataSize-1;
2021  everythingOk = false;
2022  }
2023 
2024  // pos3 = (dataSize+1)*3/4 - 1
2025  // pos3 >= 0 <==> dataSize >= 1/3
2026  // pos3 < (dataSize-1) <==> dataSize > 3
2027  unsigned int pos3 = (unsigned int) ( (((double) dataSize) + 1.)*3./4. - 1. );
2028  if (pos3 > (dataSize-1)) {
2029  pos3 = 0;
2030  everythingOk = false;
2031  }
2032  unsigned int pos3inc = pos3+1;
2033  if (pos3inc > (dataSize-1)) {
2034  pos3inc = dataSize-1;
2035  everythingOk = false;
2036  }
2037 
2038  double fraction1 = (((double) dataSize) + 1.)*1./4. - 1. - ((double) pos1);
2039  if (fraction1 < 0.) {
2040  fraction1 = 0.;
2041  everythingOk = false;
2042  }
2043  double fraction3 = (((double) dataSize) + 1.)*3./4. - 1. - ((double) pos3);
2044  if (fraction3 < 0.) {
2045  fraction3 = 0.;
2046  everythingOk = false;
2047  }
2048 
2049  if (everythingOk == false) {
2050  std::cerr << "In ScalarSequence<T>::subInterQuantileRange()"
2051  << ", worldRank = " << m_env.worldRank()
2052  << ": at least one adjustment was necessary"
2053  << std::endl;
2054  }
2055 
2056  //if (m_env.subDisplayFile()) {
2057  // *m_env.subDisplayFile() << "In ScalarSequence::subInterQuantileRange()"
2058  // << ", initialPos = " << initialPos
2059  // << ", this->subSequenceSize() = " << this->subSequenceSize()
2060  // << ", dataSize = " << dataSize
2061  // << ", sortedSequence.size() = " << sortedSequence.size()
2062  // << ", pos1 = " << pos1
2063  // << ", pos3 = " << pos3
2064  // << std::endl;
2065  //}
2066 
2067  T value1 = (1.-fraction1) * sortedSequence[pos1] + fraction1 * sortedSequence[pos1inc];
2068  T value3 = (1.-fraction3) * sortedSequence[pos3] + fraction3 * sortedSequence[pos3inc];
2069  T iqrValue = value3 - value1;
2070 
2071  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2072  *m_env.subDisplayFile() << "In ScalarSequence<T>::subInterQuantileRange()"
2073  << ": iqrValue = " << iqrValue
2074  << ", dataSize = " << dataSize
2075  << ", pos1 = " << pos1
2076  << ", pos3 = " << pos3
2077  << ", value1 = " << value1
2078  << ", value3 = " << value3
2079  << std::endl;
2080 
2081  // Save data only once into a separate file
2082  //std::ofstream* ofsvar = new std::ofstream(("sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
2083  //if ((ofsvar == NULL ) ||
2084  // (ofsvar->is_open() == false)) {
2085  // delete ofsvar;
2086  // ofsvar = new std::ofstream(("sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
2087 
2088  // *ofsvar << "var_sort_sub" << m_env.subIdString() << " = zeros(" << 1
2089  // << "," << dataSize
2090  // << ");"
2091  // << std::endl;
2092  // for (unsigned int j = 0; j < dataSize; ++j) {
2093  // *ofsvar << "var_sort_sub" << m_env.subIdString() << "(" << 1
2094  // << "," << j+1
2095  // << ") = " << sortedSequence[j]
2096  // << ";"
2097  // << std::endl;
2098  // }
2099  //}
2100  //delete ofsvar;
2101  }
2102 
2103  return iqrValue;
2104 }
int worldRank() const
Returns the process world rank.
Definition: Environment.C:262
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:75
ScalarSequence(const BaseEnvironment &env, unsigned int subSequenceSize, const std::string &name)
Default constructor.
unsigned int displayVerbosity() const
Definition: Environment.C:449
void subSort()
Sorts the sequence of scalars in the private attribute m_seq.
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
template<class T >
const T & QUESO::ScalarSequence< T >::subMaxPlain ( ) const

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

Definition at line 369 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().

370 {
371  if (m_subMaxPlain == NULL) {
372  if (m_subMinPlain == NULL) m_subMinPlain = new T(0.);
373  m_subMaxPlain = new T(0.);
375  }
376 
377  return *m_subMaxPlain;
378 }
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 838 of file ScalarSequence.C.

References queso_require_msg.

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

841 {
842  if (this->subSequenceSize() == 0) return 0.;
843 
844  bool bRC = ((initialPos < this->subSequenceSize()) &&
845  (0 < numPos ) &&
846  ((initialPos+numPos) <= this->subSequenceSize()));
847  if (bRC == false) {
848  std::cerr << "In ScalarSequence<T>::subMeanExtra()"
849  << ": ERROR at fullRank " << m_env.fullRank()
850  << ", initialPos = " << initialPos
851  << ", numPos = " << numPos
852  << ", this->subSequenceSize() = " << this->subSequenceSize()
853  << std::endl;
854  if (m_env.subDisplayFile()) {
855  *m_env.subDisplayFile() << "In ScalarSequence<T>::subMeanExtra()"
856  << ": ERROR at fullRank " << m_env.fullRank()
857  << ", initialPos = " << initialPos
858  << ", numPos = " << numPos
859  << ", this->subSequenceSize() = " << this->subSequenceSize()
860  << std::endl;
861  }
862  }
863  queso_require_msg(bRC, "invalid input data");
864 
865  unsigned int finalPosPlus1 = initialPos + numPos;
866  T tmpSum = 0.;
867  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
868  tmpSum += m_seq[j];
869  }
870 
871  return tmpSum/(T) numPos;
872 }
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
int fullRank() const
Returns the process full rank.
Definition: Environment.C:268
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
template<class T >
const T & QUESO::ScalarSequence< T >::subMeanPlain ( ) const

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

Definition at line 395 of file ScalarSequence.C.

396 {
397  if (m_subMeanPlain == NULL) {
398  m_subMeanPlain = new T(0.);
400  }
401 
402  return *m_subMeanPlain;
403 }
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 954 of file ScalarSequence.C.

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

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

957 {
958  if (this->subSequenceSize() == 0) return 0.;
959 
960  bool bRC = ((initialPos < this->subSequenceSize()) &&
961  (0 < numPos ) &&
962  ((initialPos+numPos) <= this->subSequenceSize()));
963  if (bRC == false) {
964  std::cerr << "In ScalarSequence<T>::subMedianExtra()"
965  << ": ERROR at fullRank " << m_env.fullRank()
966  << ", initialPos = " << initialPos
967  << ", numPos = " << numPos
968  << ", this->subSequenceSize() = " << this->subSequenceSize()
969  << std::endl;
970  if (m_env.subDisplayFile()) {
971  *m_env.subDisplayFile() << "In ScalarSequence<T>::subMedianExtra()"
972  << ": ERROR at fullRank " << m_env.fullRank()
973  << ", initialPos = " << initialPos
974  << ", numPos = " << numPos
975  << ", this->subSequenceSize() = " << this->subSequenceSize()
976  << std::endl;
977  }
978  }
979  queso_require_msg(bRC, "invalid input data");
980 
981  ScalarSequence sortedSequence(m_env,0,"");
982  sortedSequence.resizeSequence(numPos);
983  this->extractScalarSeq(initialPos,
984  1,
985  numPos,
986  sortedSequence);
987  sortedSequence.subSort();
988 
989  unsigned int tmpPos = (unsigned int) (0.5 * (double) numPos);
990  T resultValue = sortedSequence[tmpPos];
991 
992  return resultValue;
993 }
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
ScalarSequence(const BaseEnvironment &env, unsigned int subSequenceSize, const std::string &name)
Default constructor.
int fullRank() const
Returns the process full rank.
Definition: Environment.C:268
void extractScalarSeq(unsigned int initialPos, unsigned int spacing, unsigned int numPos, ScalarSequence< T > &scalarSeq) const
Extracts a sequence of scalars.
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
template<class T >
const T & QUESO::ScalarSequence< T >::subMedianPlain ( ) const

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

Definition at line 419 of file ScalarSequence.C.

420 {
421  if (m_subMedianPlain == NULL) {
422  m_subMedianPlain = new T(0.);
424  }
425 
426  return *m_subMedianPlain;
427 }
T subMedianExtra(unsigned int initialPos, unsigned int numPos) const
Finds the median value of the sub-sequence, considering numPos positions starting at position initial...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1501 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().

1506 {
1507  queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid input");
1508 
1510  std::advance(pos1,initialPos);
1511 
1513  std::advance(pos2,initialPos+numPos);
1514 
1515  if ((initialPos+numPos) == this->subSequenceSize()) {
1516  queso_require_msg(!(pos2 != m_seq.end()), "invalid state");
1517  }
1518 
1520  pos = std::min_element(pos1, pos2);
1521  minValue = *pos;
1522  pos = std::max_element(pos1, pos2);
1523  maxValue = *pos;
1524 
1525  return;
1526 }
std::vector< T >::const_iterator seqScalarPositionConstIteratorTypedef
#define queso_require_less_equal_msg(expr1, expr2, msg)
Definition: asserts.h:77
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
template<class T >
const T & QUESO::ScalarSequence< T >::subMinPlain ( ) const

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

Definition at line 343 of file ScalarSequence.C.

344 {
345  if (m_subMinPlain == NULL) {
346  m_subMinPlain = new T(0.);
347  if (m_subMaxPlain == NULL) m_subMaxPlain = new T(0.);
349  }
350 
351  return *m_subMinPlain;
352 }
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1218 of file ScalarSequence.C.

References queso_require_msg.

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

1222 {
1223  if (this->subSequenceSize() == 0) return 0.;
1224 
1225  bool bRC = ((initialPos < this->subSequenceSize()) &&
1226  (0 < numPos ) &&
1227  ((initialPos+numPos) <= this->subSequenceSize()));
1228  queso_require_msg(bRC, "invalid input data");
1229 
1230  unsigned int finalPosPlus1 = initialPos + numPos;
1231  T diff;
1232  T popValue = 0.;
1233  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1234  diff = m_seq[j] - meanValue;
1235  popValue += diff*diff;
1236  }
1237 
1238  popValue /= (T) numPos;
1239 
1240  return popValue;
1241 }
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2509 of file ScalarSequence.C.

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

2512 {
2513  queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
2514 
2515  T subMaxValue = subCorrespondingScalarValues.subMaxPlain();
2516  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
2517 
2518  unsigned int subNumPos = 0;
2519  for (unsigned int i = 0; i < iMax; ++i) {
2520  if (subCorrespondingScalarValues[i] == subMaxValue) {
2521  subNumPos++;
2522  }
2523  }
2524 
2525  subPositionsOfMaximum.resizeSequence(subNumPos);
2526  unsigned int j = 0;
2527  for (unsigned int i = 0; i < iMax; ++i) {
2528  if (subCorrespondingScalarValues[i] == subMaxValue) {
2529  subPositionsOfMaximum[j] = (*this)[i];
2530  j++;
2531  }
2532  }
2533 
2534  return subMaxValue;
2535 }
T subPositionsOfMaximum(const ScalarSequence< T > &subCorrespondingScalarValues, ScalarSequence< T > &subPositionsOfMaximum)
Finds the positions where the maximum element occurs in the sub-sequence.
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1130 of file ScalarSequence.C.

References queso_require_msg.

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

1134 {
1135  if (this->subSequenceSize() == 0) return 0.;
1136 
1137  bool bRC = ((initialPos < this->subSequenceSize()) &&
1138  (0 < numPos ) &&
1139  ((initialPos+numPos) <= this->subSequenceSize()));
1140  queso_require_msg(bRC, "invalid input data");
1141 
1142  unsigned int finalPosPlus1 = initialPos + numPos;
1143  T diff;
1144  T stdValue = 0.;
1145  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1146  diff = m_seq[j] - meanValue;
1147  stdValue += diff*diff;
1148  }
1149 
1150  stdValue /= (((T) numPos) - 1.);
1151  stdValue = sqrt(stdValue);
1152 
1153  return stdValue;
1154 }
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1044 of file ScalarSequence.C.

References queso_require_msg.

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

1048 {
1049  if (this->subSequenceSize() == 0) return 0.;
1050 
1051  bool bRC = ((initialPos < this->subSequenceSize()) &&
1052  (0 < numPos ) &&
1053  ((initialPos+numPos) <= this->subSequenceSize()));
1054  queso_require_msg(bRC, "invalid input data");
1055 
1056  unsigned int finalPosPlus1 = initialPos + numPos;
1057  T diff;
1058  T samValue = 0.;
1059  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1060  diff = m_seq[j] - meanValue;
1061  samValue += diff*diff;
1062  }
1063 
1064  samValue /= (((T) numPos) - 1.);
1065 
1066  return samValue;
1067 }
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 443 of file ScalarSequence.C.

444 {
445  if (m_subSampleVariancePlain == NULL) {
446  m_subSampleVariancePlain = new T(0.);
448  }
449 
450  return *m_subSampleVariancePlain;
451 }
const T & subMeanPlain() const
Finds the mean value of the sub-sequence of scalars.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 init...
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 2199 of file ScalarSequence.C.

References queso_require_msg.

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

2203 {
2204  bool bRC = (initialPos < this->subSequenceSize());
2205  queso_require_msg(bRC, "invalid input data");
2206 
2207  unsigned int dataSize = this->subSequenceSize() - initialPos;
2208 
2209  T meanValue = this->subMeanExtra(initialPos,
2210  dataSize);
2211 
2212  T samValue = this->subSampleVarianceExtra(initialPos,
2213  dataSize,
2214  meanValue);
2215 
2216  T scaleValue;
2217  if (iqrValue <= 0.) {
2218  scaleValue = 1.06*std::sqrt(samValue)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
2219  }
2220  else {
2221  scaleValue = 1.06*std::min(std::sqrt(samValue),iqrValue/1.34)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
2222  }
2223 
2224  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2225  *m_env.subDisplayFile() << "In ScalarSequence<T>::subScaleForKde()"
2226  << ": iqrValue = " << iqrValue
2227  << ", meanValue = " << meanValue
2228  << ", samValue = " << samValue
2229  << ", dataSize = " << dataSize
2230  << ", scaleValue = " << scaleValue
2231  << std::endl;
2232  }
2233 
2234  return scaleValue;
2235 }
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
unsigned int displayVerbosity() const
Definition: Environment.C:449
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 init...
template<class T >
unsigned int QUESO::ScalarSequence< T >::subSequenceSize ( ) const
template<class T>
void QUESO::ScalarSequence< T >::subSort ( unsigned int  initialPos,
ScalarSequence< T > &  sortedSequence 
) const

Sorts the sub-sequence of scalars.

Definition at line 1874 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().

1877 {
1878  unsigned int numPos = this->subSequenceSize() - initialPos;
1879  sortedSequence.resizeSequence(numPos);
1880  this->extractScalarSeq(initialPos,
1881  1,
1882  numPos,
1883  sortedSequence);
1884  sortedSequence.subSort();
1885 
1886  return;
1887 }
void extractScalarSeq(unsigned int initialPos, unsigned int spacing, unsigned int numPos, ScalarSequence< T > &scalarSeq) const
Extracts a sequence of scalars.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
template<class T>
void QUESO::ScalarSequence< T >::subSort ( )
private

Sorts the sequence of scalars in the private attribute m_seq.

Definition at line 3255 of file ScalarSequence.C.

3256 {
3257  std::sort(m_seq.begin(), m_seq.end());
3258  return;
3259 }
std::vector< T > m_seq
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 571 of file ScalarSequence.C.

576 {
577  T tmpMinValue;
578  T tmpMaxValue;
579  std::vector<T> centers(numEvaluationPoints,0.);
580  std::vector<unsigned int> bins (numEvaluationPoints,0);
581 
582  subMinMaxExtra(0, // initialPos
583  this->subSequenceSize(),
584  tmpMinValue,
585  tmpMaxValue);
586  subHistogram(0, // initialPos,
587  tmpMinValue,
588  tmpMaxValue,
589  centers,
590  bins);
591 
592  minDomainValue = centers[0];
593  maxDomainValue = centers[centers.size()-1];
594 
595  unsigned int sumOfBins = 0;
596  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
597  sumOfBins += bins[i];
598  }
599 
600  cdfValues.clear();
601  cdfValues.resize(numEvaluationPoints);
602  unsigned int partialSum = 0;
603  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
604  partialSum += bins[i];
605  cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
606  }
607 
608  return;
609 }
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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.
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 749 of file ScalarSequence.C.

753 {
754  T tmpMinValue;
755  T tmpMaxValue;
756  std::vector<unsigned int> bins(numEvaluationPoints,0);
757  gridValues.resize (numEvaluationPoints,0.);
758  cdfValues.resize (numEvaluationPoints,0.);
759 
760  subMinMaxExtra(0, // initialPos
761  this->subSequenceSize(),
762  tmpMinValue,
763  tmpMaxValue);
764 
765  if (tmpMinValue == tmpMaxValue) {
766  if (tmpMinValue < -1.e-12) {
767  tmpMinValue += tmpMinValue*(1.e-8);
768  }
769  else if (tmpMinValue > 1.e-12) {
770  tmpMinValue -= tmpMinValue*(1.e-8);
771  }
772  else {
773  tmpMinValue = 1.e-12;
774  }
775  }
776 
777  subWeightHistogram(0, // initialPos,
778  tmpMinValue,
779  tmpMaxValue,
780  gridValues,
781  bins);
782 
783  unsigned int sumOfBins = 0;
784  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
785  sumOfBins += bins[i];
786  }
787 
788  cdfValues.clear();
789  cdfValues.resize(numEvaluationPoints);
790  unsigned int partialSum = 0;
791  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
792  partialSum += bins[i];
793  cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
794  }
795 
796  return;
797 }
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.
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 801 of file ScalarSequence.C.

805 {
806  T tmpMinValue;
807  T tmpMaxValue;
808  std::vector<unsigned int> bins(numEvaluationPoints,0);
809 
810  subMinMaxExtra(0, // initialPos
811  this->subSequenceSize(),
812  tmpMinValue,
813  tmpMaxValue);
814  subWeightHistogram(0, // initialPos,
815  tmpMinValue,
816  tmpMaxValue,
817  gridValues,
818  bins);
819 
820  unsigned int sumOfBins = 0;
821  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
822  sumOfBins += bins[i];
823  }
824 
825  cdfValues.clear();
826  cdfValues.resize(numEvaluationPoints);
827  unsigned int partialSum = 0;
828  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
829  partialSum += bins[i];
830  cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
831  }
832 
833  return;
834 }
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.
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1780 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

1786 {
1787  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1788 
1789  for (unsigned int j = 0; j < bins.size(); ++j) {
1790  bins[j] = 0;
1791  }
1792 
1793  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1794  double minCenter = minHorizontalValue - horizontalDelta/2.;
1795  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1796  gridValues = new UniformOneDGrid<T>(m_env,
1797  "",
1798  bins.size(),
1799  minCenter,
1800  maxCenter);
1801 
1802  unsigned int dataSize = this->subSequenceSize();
1803  for (unsigned int j = 0; j < dataSize; ++j) {
1804  double value = m_seq[j];
1805  if (value < minHorizontalValue) {
1806  bins[0] += value;
1807  }
1808  else if (value >= maxHorizontalValue) {
1809  bins[bins.size()-1] += value;
1810  }
1811  else {
1812  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1813  bins[index] += value;
1814  }
1815  }
1816 
1817  return;
1818 }
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
std::vector< T > m_seq
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1822 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

1828 {
1829  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1830 
1831  for (unsigned int j = 0; j < bins.size(); ++j) {
1832  bins[j] = 0;
1833  }
1834 
1835  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1836  double minCenter = minHorizontalValue - horizontalDelta/2.;
1837  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1838  UniformOneDGrid<T> tmpGrid(m_env,
1839  "",
1840  bins.size(),
1841  minCenter,
1842  maxCenter);
1843  gridValues.clear();
1844  gridValues.resize(tmpGrid.size(),0.);
1845  for (unsigned int i = 0; i < tmpGrid.size(); ++i) {
1846  gridValues[i] = tmpGrid[i];
1847  }
1848 
1849  unsigned int dataSize = this->subSequenceSize();
1850  for (unsigned int j = 0; j < dataSize; ++j) {
1851  double value = m_seq[j];
1852  if (value < minHorizontalValue) {
1853  bins[0] += value;
1854  }
1855  else if (value >= maxHorizontalValue) {
1856  bins[bins.size()-1] += value;
1857  }
1858  else {
1859  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1860  bins[index] += value;
1861  }
1862  }
1863 
1864  //std::cout << "In ScalarSequence<T>::subWeightHistogram():" << std::endl;
1865  //for (unsigned int j = 0; j < bins.size(); ++j) {
1866  // std::cout << "bins[" << j << "] = " << bins[j] << std::endl;
1867  //}
1868 
1869  return;
1870 }
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
std::vector< T > m_seq
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2569 of file ScalarSequence.C.

References queso_require_greater_equal_msg, UQ_FILE_EXTENSION_FOR_HDF_FORMAT, UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, and UQ_FILE_EXTENSION_FOR_TXT_FORMAT.

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

2575 {
2576  queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");
2577 
2578  FilePtrSetStruct filePtrSet;
2579  if (m_env.openOutputFile(fileName,
2580  fileType,
2581  allowedSubEnvIds,
2582  false, // A 'true' causes problems when the user chooses (via options
2583  // in the input file) to use just one file for all outputs.
2584  filePtrSet)) {
2585 
2586  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT ||
2587  fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) {
2588  this->subWriteContents(initialPos,
2589  numPos,
2590  *filePtrSet.ofsVar,
2591  fileType);
2592  }
2593 #ifdef QUESO_HAS_HDF5
2594  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2595 
2596  // Create dataspace
2597  hsize_t dims[1] = { this->subSequenceSize() };
2598  hid_t dataspace_id = H5Screate_simple(1, dims, dims);
2600  dataspace_id,
2601  0,
2602  "error create dataspace with id: " << dataspace_id);
2603 
2604  // Create dataset
2605  hid_t dataset_id = H5Dcreate(filePtrSet.h5Var,
2606  "data",
2607  H5T_IEEE_F64LE,
2608  dataspace_id,
2609  H5P_DEFAULT,
2610  H5P_DEFAULT,
2611  H5P_DEFAULT);
2612 
2614  dataset_id,
2615  0,
2616  "error creating dataset with id: " << dataset_id);
2617 
2618  // Write the dataset
2619  herr_t status = H5Dwrite(
2620  dataset_id,
2621  H5T_NATIVE_DOUBLE, // The type in memory
2622  H5S_ALL, // The dataspace in memory
2623  dataspace_id, // The file dataspace
2624  H5P_DEFAULT, // Xfer property list
2625  &m_seq[0]);
2626 
2628  status,
2629  0,
2630  "error writing dataset to file with id: " << filePtrSet.h5Var);
2631 
2632  // Clean up
2633  H5Dclose(dataset_id);
2634  H5Sclose(dataspace_id);
2635  }
2636 #endif
2637 
2638  m_env.closeFile(filePtrSet,fileType);
2639  }
2640 }
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
int subRank() const
Access function for sub-rank.
Definition: Environment.C:287
#define UQ_FILE_EXTENSION_FOR_HDF_FORMAT
Definition: Defines.h:104
#define UQ_FILE_EXTENSION_FOR_TXT_FORMAT
Definition: Defines.h:103
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.
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
Definition: Environment.C:1083
bool openOutputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, bool writeOver, FilePtrSetStruct &filePtrSet) const
Opens an output file for each sub-environment that was chosen to send data to the file...
Definition: Environment.C:520
std::vector< T > m_seq
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:102
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2644 of file ScalarSequence.C.

References UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, and UQ_FILE_EXTENSION_FOR_TXT_FORMAT.

2649 {
2650  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2651  this->writeSubMatlabHeader(ofs, this->subSequenceSize());
2652  }
2653  else if (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) {
2654  this->writeTxtHeader(ofs, this->subSequenceSize());
2655  }
2656 
2657  unsigned int chainSize = this->subSequenceSize();
2658  for (unsigned int j = 0; j < chainSize; ++j) {
2659  ofs << m_seq[j]
2660  << std::endl;
2661  }
2662 
2663  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2664  ofs << "];\n";
2665  }
2666 
2667  return;
2668 }
void writeSubMatlabHeader(std::ofstream &ofs, double sequenceSize) const
Helper function to write header info for matlab files from one chain.
#define UQ_FILE_EXTENSION_FOR_TXT_FORMAT
Definition: Defines.h:103
std::vector< T > m_seq
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:102
void writeTxtHeader(std::ofstream &ofs, double sequenceSize) const
Helper function to write txt info for matlab files.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2340 of file ScalarSequence.C.

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

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

2346 {
2347  if (m_env.numSubEnvironments() == 1) {
2348  return this->subGaussian1dKde(initialPos,
2349  unifiedScaleValue,
2350  unifiedEvaluationPositions,
2351  unifiedDensityValues);
2352  }
2353 
2354  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
2355 
2356  if (useOnlyInter0Comm) {
2357  if (m_env.inter0Rank() >= 0) {
2358  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2359  (0 < unifiedEvaluationPositions.size()) &&
2360  (unifiedEvaluationPositions.size() == unifiedDensityValues.size() ));
2361  queso_require_msg(bRC, "invalid input data");
2362 
2363  unsigned int localDataSize = this->subSequenceSize() - initialPos;
2364  unsigned int unifiedDataSize = 0;
2365  m_env.inter0Comm().template Allreduce<unsigned int>(&localDataSize, &unifiedDataSize, (int) 1, RawValue_MPI_SUM,
2366  "ScalarSequence<T>::unifiedGaussian1dKde()",
2367  "failed MPI.Allreduce() for data size");
2368 
2369  unsigned int numEvals = unifiedEvaluationPositions.size();
2370 
2371  std::vector<double> densityValues(numEvals,0.);
2372  double unifiedScaleInv = 1./unifiedScaleValue;
2373  for (unsigned int j = 0; j < numEvals; ++j) {
2374  double x = unifiedEvaluationPositions[j];
2375  double value = 0.;
2376  for (unsigned int k = 0; k < localDataSize; ++k) {
2377  double xk = m_seq[initialPos+k];
2378  value += MiscGaussianDensity((x-xk)*unifiedScaleInv,0.,1.);
2379  }
2380  densityValues[j] = value;
2381  }
2382 
2383  for (unsigned int j = 0; j < numEvals; ++j) {
2384  unifiedDensityValues[j] = 0.;
2385  }
2386  m_env.inter0Comm().template Allreduce<double>(&densityValues[0], &unifiedDensityValues[0], (int) numEvals, RawValue_MPI_SUM,
2387  "ScalarSequence<T>::unifiedGaussian1dKde()",
2388  "failed MPI.Allreduce() for density values");
2389 
2390  for (unsigned int j = 0; j < numEvals; ++j) {
2391  unifiedDensityValues[j] *= unifiedScaleInv/((double) unifiedDataSize);
2392  }
2393 
2394  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2395  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedGaussian1dKde()"
2396  << ": unifiedDensityValues[0] = " << unifiedDensityValues[0]
2397  << ", unifiedDensityValues[" << unifiedDensityValues.size()-1 << "] = " << unifiedDensityValues[unifiedDensityValues.size()-1]
2398  << std::endl;
2399  }
2400  }
2401  else {
2402  // Node not in the 'inter0' communicator
2403  this->subGaussian1dKde(initialPos,
2404  unifiedScaleValue,
2405  unifiedEvaluationPositions,
2406  unifiedDensityValues);
2407  }
2408  }
2409  else {
2410  queso_error_msg("parallel vectors not supported yet");
2411  }
2412 
2413  //m_env.fullComm().Barrier();
2414 
2415  return;
2416 }
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
double MiscGaussianDensity(double x, double mu, double sigma)
int k
Definition: ann_sample.cpp:53
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.
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
unsigned int displayVerbosity() const
Definition: Environment.C:449
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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 1648 of file ScalarSequence.C.

References queso_error_msg, queso_require_equal_to_msg, queso_require_greater_equal_msg, and RawValue_MPI_SUM.

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

1655 {
1656  if (m_env.numSubEnvironments() == 1) {
1657  return this->subHistogram(initialPos,
1658  unifiedMinHorizontalValue,
1659  unifiedMaxHorizontalValue,
1660  unifiedCenters,
1661  unifiedBins);
1662  }
1663 
1664  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
1665 
1666  if (useOnlyInter0Comm) {
1667  if (m_env.inter0Rank() >= 0) {
1668  queso_require_equal_to_msg(unifiedCenters.size(), unifiedBins.size(), "vectors 'unifiedCenters' and 'unifiedBins' have different sizes");
1669 
1670  queso_require_greater_equal_msg(unifiedBins.size(), 3, "number of 'unifiedBins' is too small: should be at least 3");
1671 
1672  for (unsigned int j = 0; j < unifiedBins.size(); ++j) {
1673  unifiedCenters[j] = 0.;
1674  unifiedBins[j] = 0;
1675  }
1676 
1677  double unifiedHorizontalDelta = (unifiedMaxHorizontalValue - unifiedMinHorizontalValue)/(((double) unifiedBins.size()) - 2.); // IMPORTANT: -2
1678 
1679  double unifiedMinCenter = unifiedMinHorizontalValue - unifiedHorizontalDelta/2.;
1680  double unifiedMaxCenter = unifiedMaxHorizontalValue + unifiedHorizontalDelta/2.;
1681  for (unsigned int j = 0; j < unifiedCenters.size(); ++j) {
1682  double factor = ((double) j)/(((double) unifiedCenters.size()) - 1.);
1683  unifiedCenters[j] = (1. - factor) * unifiedMinCenter + factor * unifiedMaxCenter;
1684  }
1685 
1686  std::vector<unsigned int> localBins(unifiedBins.size(),0);
1687  unsigned int dataSize = this->subSequenceSize();
1688  for (unsigned int j = 0; j < dataSize; ++j) {
1689  double value = m_seq[j];
1690  if (value < unifiedMinHorizontalValue) {
1691  localBins[0]++;
1692  }
1693  else if (value >= unifiedMaxHorizontalValue) {
1694  localBins[localBins.size()-1]++;
1695  }
1696  else {
1697  unsigned int index = 1 + (unsigned int) ((value - unifiedMinHorizontalValue)/unifiedHorizontalDelta);
1698  localBins[index]++;
1699  }
1700  }
1701 
1702  m_env.inter0Comm().template Allreduce<unsigned int>(&localBins[0], &unifiedBins[0], (int) localBins.size(), RawValue_MPI_SUM,
1703  "ScalarSequence<T>::unifiedHistogram()",
1704  "failed MPI.Allreduce() for bins");
1705 
1706  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1707  for (unsigned int i = 0; i < unifiedCenters.size(); ++i) {
1708  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedHistogram()"
1709  << ": i = " << i
1710  << ", unifiedMinHorizontalValue = " << unifiedMinHorizontalValue
1711  << ", unifiedMaxHorizontalValue = " << unifiedMaxHorizontalValue
1712  << ", unifiedCenters = " << unifiedCenters[i]
1713  << ", unifiedBins = " << unifiedBins[i]
1714  << std::endl;
1715  }
1716  }
1717  }
1718  else {
1719  // Node not in the 'inter0' communicator
1720  this->subHistogram(initialPos,
1721  unifiedMinHorizontalValue,
1722  unifiedMaxHorizontalValue,
1723  unifiedCenters,
1724  unifiedBins);
1725  }
1726  }
1727  else {
1728  queso_error_msg("parallel vectors not supported yet");
1729  }
1730 
1731  //m_env.fullComm().Barrier();
1732 
1733  return;
1734 }
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
unsigned int displayVerbosity() const
Definition: Environment.C:449
std::vector< T > m_seq
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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.
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 2108 of file ScalarSequence.C.

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

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

2111 {
2112  T unifiedIqrValue = 0.;
2113 
2114  if (m_env.numSubEnvironments() == 1) {
2115  return this->subInterQuantileRange(initialPos);
2116  }
2117 
2118  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
2119 
2120  if (useOnlyInter0Comm) {
2121  if (m_env.inter0Rank() >= 0) {
2122  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedInterQuantileRange(), beginning logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
2123 
2124  ScalarSequence unifiedSortedSequence(m_env,0,"");
2125  this->unifiedSort(useOnlyInter0Comm,
2126  initialPos,
2127  unifiedSortedSequence);
2128  unsigned int unifiedDataSize = unifiedSortedSequence.subSequenceSize();
2129 
2130  unsigned int localDataSize = this->subSequenceSize() - initialPos;
2131  unsigned int sumOfLocalSizes = 0;
2132  m_env.inter0Comm().template Allreduce<unsigned int>(&localDataSize, &sumOfLocalSizes, (int) 1, RawValue_MPI_SUM,
2133  "ScalarSequence<T>::unifiedInterQuantileRange()",
2134  "failed MPI.Allreduce() for data size");
2135 
2136  queso_require_equal_to_msg(sumOfLocalSizes, unifiedDataSize, "incompatible unified sizes");
2137 
2138  unsigned int pos1 = (unsigned int) ( (((double) unifiedDataSize) + 1.)*1./4. - 1. );
2139  unsigned int pos3 = (unsigned int) ( (((double) unifiedDataSize) + 1.)*3./4. - 1. );
2140 
2141  double fraction1 = (((double) unifiedDataSize) + 1.)*1./4. - 1. - ((double) pos1);
2142  double fraction3 = (((double) unifiedDataSize) + 1.)*3./4. - 1. - ((double) pos3);
2143 
2144  T value1 = (1.-fraction1) * unifiedSortedSequence[pos1] + fraction1 * unifiedSortedSequence[pos1+1];
2145  T value3 = (1.-fraction3) * unifiedSortedSequence[pos3] + fraction3 * unifiedSortedSequence[pos3+1];
2146  unifiedIqrValue = value3 - value1;
2147 
2148  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2149  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedInterQuantileRange()"
2150  << ": unifiedIqrValue = " << unifiedIqrValue
2151  << ", localDataSize = " << localDataSize
2152  << ", unifiedDataSize = " << unifiedDataSize
2153  << ", pos1 = " << pos1
2154  << ", pos3 = " << pos3
2155  << ", value1 = " << value1
2156  << ", value3 = " << value3
2157  << std::endl;
2158 
2159  // Save data only once into a separate file
2160  //std::ofstream* ofsvar = new std::ofstream(("unif_sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
2161  //if ((ofsvar == NULL ) ||
2162  // (ofsvar->is_open() == false)) {
2163  // delete ofsvar;
2164  // ofsvar = new std::ofstream(("unif_sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
2165 
2166  // *ofsvar << "var_unif_sort_sub" << m_env.subIdString() << " = zeros(" << 1
2167  // << "," << unifiedDataSize
2168  // << ");"
2169  // << std::endl;
2170  // for (unsigned int j = 0; j < unifiedDataSize; ++j) {
2171  // *ofsvar << "var_unif_sort_sub" << m_env.subIdString() << "(" << 1
2172  // << "," << j+1
2173  // << ") = " << unifiedSortedSequence[j]
2174  // << ";"
2175  // << std::endl;
2176  // }
2177  //}
2178  //delete ofsvar;
2179  }
2180 
2181  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedInterQuantileRange(), ending logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
2182  }
2183  else {
2184  // Node not in the 'inter0' communicator
2185  unifiedIqrValue = this->subInterQuantileRange(initialPos);
2186  }
2187  }
2188  else {
2189  queso_error_msg("parallel vectors not supported yet");
2190  }
2191 
2192  //m_env.fullComm().Barrier();
2193 
2194  return unifiedIqrValue;
2195 }
void unifiedSort(bool useOnlyInter0Comm, unsigned int initialPos, ScalarSequence< T > &unifiedSortedSequence) const
Sorts the unified sequence of scalars.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
ScalarSequence(const BaseEnvironment &env, unsigned int subSequenceSize, const std::string &name)
Default constructor.
T subInterQuantileRange(unsigned int initialPos) const
Returns the interquartile range of the values in the sub-sequence.
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
unsigned int displayVerbosity() const
Definition: Environment.C:449
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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 382 of file ScalarSequence.C.

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

383 {
384  if (m_unifiedMaxPlain == NULL) {
385  if (m_unifiedMinPlain == NULL) m_unifiedMinPlain = new T(0.);
386  m_unifiedMaxPlain = new T(0.);
388  }
389 
390  return *m_unifiedMaxPlain;
391 }
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 starti...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 876 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, and RawValue_MPI_SUM.

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

880 {
881  if (m_env.numSubEnvironments() == 1) {
882  return this->subMeanExtra(initialPos,
883  numPos);
884  }
885 
886  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
887 
888  T unifiedMeanValue = 0.;
889  if (useOnlyInter0Comm) {
890  if (m_env.inter0Rank() >= 0) {
891  bool bRC = ((initialPos < this->subSequenceSize()) &&
892  (0 < numPos ) &&
893  ((initialPos+numPos) <= this->subSequenceSize()));
894  queso_require_msg(bRC, "invalid input data");
895 
896  unsigned int finalPosPlus1 = initialPos + numPos;
897  T localSum = 0.;
898  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
899  localSum += m_seq[j];
900  }
901 
902  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
903  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMeanExtra()"
904  << ": initialPos = " << initialPos
905  << ", numPos = " << numPos
906  << ", before MPI.Allreduce"
907  << std::endl;
908  }
909  //std::cout << m_env.inter0Comm().MyPID()
910  // << std::endl;
911  //sleep(1);
912  unsigned int unifiedNumPos = 0;
913  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
914  "ScalarSequence<T>::unifiedMeanExtra()",
915  "failed MPI.Allreduce() for numPos");
916 
917  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
918  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMeanExtra()"
919  << ": numPos = " << numPos
920  << ", unifiedNumPos = " << unifiedNumPos
921  << std::endl;
922  }
923 
924  m_env.inter0Comm().template Allreduce<double>(&localSum, &unifiedMeanValue, (int) 1, RawValue_MPI_SUM,
925  "ScalarSequence<T>::unifiedMeanExtra()",
926  "failed MPI.Allreduce() for sum");
927 
928  unifiedMeanValue /= ((T) unifiedNumPos);
929 
930  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
931  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMeanExtra()"
932  << ": localSum = " << localSum
933  << ", unifiedMeanValue = " << unifiedMeanValue
934  << std::endl;
935  }
936  }
937  else {
938  // Node not in the 'inter0' communicator
939  this->subMeanExtra(initialPos,
940  numPos);
941  }
942  }
943  else {
944  queso_error_msg("parallel vectors not supported yet");
945  }
946 
947  //m_env.fullComm().Barrier();
948 
949  return unifiedMeanValue;
950 }
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
unsigned int displayVerbosity() const
Definition: Environment.C:449
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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 407 of file ScalarSequence.C.

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

408 {
409  if (m_unifiedMeanPlain == NULL) {
410  m_unifiedMeanPlain = new T(0.);
411  *m_unifiedMeanPlain = unifiedMeanExtra(useOnlyInter0Comm,0,subSequenceSize());
412  }
413 
414  return *m_unifiedMeanPlain;
415 }
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...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 997 of file ScalarSequence.C.

References queso_error_msg, and queso_require_msg.

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

1001 {
1002  if (m_env.numSubEnvironments() == 1) {
1003  return this->subMedianExtra(initialPos,
1004  numPos);
1005  }
1006 
1007  // As of 07/Jul/2012, this routine does *not* require sub sequences to have equal size. Good.
1008 
1009  T unifiedMedianValue = 0.;
1010  if (useOnlyInter0Comm) {
1011  if (m_env.inter0Rank() >= 0) {
1012  bool bRC = ((initialPos < this->subSequenceSize()) &&
1013  (0 < numPos ) &&
1014  ((initialPos+numPos) <= this->subSequenceSize()));
1015  queso_require_msg(bRC, "invalid input data");
1016 
1017  ScalarSequence unifiedSortedSequence(m_env,0,"");
1018  this->unifiedSort(useOnlyInter0Comm,
1019  initialPos,
1020  unifiedSortedSequence);
1021  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1022  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMedianExtra()"
1023  << ", unifiedMedianValue = " << unifiedMedianValue
1024  << std::endl;
1025  }
1026  }
1027  else {
1028  // Node not in the 'inter0' communicator
1029  this->subMedianExtra(initialPos,
1030  numPos);
1031  }
1032  }
1033  else {
1034  queso_error_msg("parallel vectors not supported yet");
1035  }
1036 
1037  //m_env.fullComm().Barrier();
1038 
1039  return unifiedMedianValue;
1040 }
void unifiedSort(bool useOnlyInter0Comm, unsigned int initialPos, ScalarSequence< T > &unifiedSortedSequence) const
Sorts the unified sequence of scalars.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
T subMedianExtra(unsigned int initialPos, unsigned int numPos) const
Finds the median value of the sub-sequence, considering numPos positions starting at position initial...
ScalarSequence(const BaseEnvironment &env, unsigned int subSequenceSize, const std::string &name)
Default constructor.
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
unsigned int displayVerbosity() const
Definition: Environment.C:449
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 431 of file ScalarSequence.C.

432 {
433  if (m_unifiedMedianPlain == NULL) {
434  m_unifiedMedianPlain = new T(0.);
435  *m_unifiedMedianPlain = unifiedMedianExtra(useOnlyInter0Comm,0,subSequenceSize());
436  }
437 
438  return *m_unifiedMedianPlain;
439 }
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 ini...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1530 of file ScalarSequence.C.

References queso_error_msg, RawValue_MPI_MAX, and RawValue_MPI_MIN.

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

1536 {
1537  if (m_env.numSubEnvironments() == 1) {
1538  return this->subMinMaxExtra(initialPos,
1539  numPos,
1540  unifiedMinValue,
1541  unifiedMaxValue);
1542  }
1543 
1544  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1545 
1546  if (useOnlyInter0Comm) {
1547  if (m_env.inter0Rank() >= 0) {
1548  // Find local min and max
1549  T minValue;
1550  T maxValue;
1551  this->subMinMaxExtra(initialPos,
1552  numPos,
1553  minValue,
1554  maxValue);
1555 
1556  // Get overall min
1557  std::vector<double> sendBuf(1,0.);
1558  for (unsigned int i = 0; i < sendBuf.size(); ++i) {
1559  sendBuf[i] = minValue;
1560  }
1561  m_env.inter0Comm().template Allreduce<double>(&sendBuf[0], &unifiedMinValue, (int) sendBuf.size(), RawValue_MPI_MIN,
1562  "ScalarSequence<T>::unifiedMinMaxExtra()",
1563  "failed MPI.Allreduce() for min");
1564 
1565  // Get overall max
1566  for (unsigned int i = 0; i < sendBuf.size(); ++i) {
1567  sendBuf[i] = maxValue;
1568  }
1569  m_env.inter0Comm().template Allreduce<double>(&sendBuf[0], &unifiedMaxValue, (int) sendBuf.size(), RawValue_MPI_MAX,
1570  "ScalarSequence<T>::unifiedMinMaxExtra()",
1571  "failed MPI.Allreduce() for max");
1572 
1573  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1574  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMinMaxExtra()"
1575  << ": localMinValue = " << minValue
1576  << ", localMaxValue = " << maxValue
1577  << ", unifiedMinValue = " << unifiedMinValue
1578  << ", unifiedMaxValue = " << unifiedMaxValue
1579  << std::endl;
1580  }
1581  }
1582  else {
1583  // Node not in the 'inter0' communicator
1584  this->subMinMaxExtra(initialPos,
1585  numPos,
1586  unifiedMinValue,
1587  unifiedMaxValue);
1588  }
1589  }
1590  else {
1591  queso_error_msg("parallel vectors not supported yet");
1592  }
1593 
1594  //m_env.fullComm().Barrier();
1595 
1596  return;
1597 }
#define RawValue_MPI_MIN
Definition: MpiComm.h:69
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
#define RawValue_MPI_MAX
Definition: MpiComm.h:70
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
unsigned int displayVerbosity() const
Definition: Environment.C:449
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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 356 of file ScalarSequence.C.

357 {
358  if (m_unifiedMinPlain == NULL) {
359  m_unifiedMinPlain = new T(0.);
360  if (m_unifiedMaxPlain == NULL) m_unifiedMaxPlain = new T(0.);
362  }
363 
364  return *m_unifiedMinPlain;
365 }
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 starti...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 1245 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, and RawValue_MPI_SUM.

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

1250 {
1251  if (m_env.numSubEnvironments() == 1) {
1252  return this->subPopulationVariance(initialPos,
1253  numPos,
1254  unifiedMeanValue);
1255  }
1256 
1257  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1258 
1259  T unifiedPopValue = 0.;
1260  if (useOnlyInter0Comm) {
1261  if (m_env.inter0Rank() >= 0) {
1262  bool bRC = ((initialPos < this->subSequenceSize()) &&
1263  (0 < numPos ) &&
1264  ((initialPos+numPos) <= this->subSequenceSize()));
1265  queso_require_msg(bRC, "invalid input data");
1266 
1267  unsigned int finalPosPlus1 = initialPos + numPos;
1268  T diff;
1269  T localPopValue = 0.;
1270  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1271  diff = m_seq[j] - unifiedMeanValue;
1272  localPopValue += diff*diff;
1273  }
1274 
1275  unsigned int unifiedNumPos = 0;
1276  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
1277  "ScalarSequence<T>::unifiedPopulationVariance()",
1278  "failed MPI.Allreduce() for numPos");
1279 
1280  m_env.inter0Comm().template Allreduce<double>(&localPopValue, &unifiedPopValue, (int) 1, RawValue_MPI_SUM,
1281  "ScalarSequence<T>::unifiedPopulationVariance()",
1282  "failed MPI.Allreduce() for popValue");
1283 
1284  unifiedPopValue /= ((T) unifiedNumPos);
1285  }
1286  else {
1287  // Node not in the 'inter0' communicator
1288  this->subPopulationVariance(initialPos,
1289  numPos,
1290  unifiedMeanValue);
1291  }
1292  }
1293  else {
1294  queso_error_msg("parallel vectors not supported yet");
1295  }
1296 
1297  //m_env.fullComm().Barrier();
1298 
1299  return unifiedPopValue;
1300 }
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
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 ...
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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 2539 of file ScalarSequence.C.

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

2542 {
2543  queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
2544 
2545  T maxValue = subCorrespondingScalarValues.subMaxPlain();
2546  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
2547 
2548  unsigned int numPos = 0;
2549  for (unsigned int i = 0; i < iMax; ++i) {
2550  if (subCorrespondingScalarValues[i] == maxValue) {
2551  numPos++;
2552  }
2553  }
2554 
2555  unifiedPositionsOfMaximum.resizeSequence(numPos);
2556  unsigned int j = 0;
2557  for (unsigned int i = 0; i < iMax; ++i) {
2558  if (subCorrespondingScalarValues[i] == maxValue) {
2559  unifiedPositionsOfMaximum[j] = (*this)[i];
2560  j++;
2561  }
2562  }
2563 
2564  return maxValue;
2565 }
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
T unifiedPositionsOfMaximum(const ScalarSequence< T > &subCorrespondingScalarValues, ScalarSequence< T > &unifiedPositionsOfMaximum)
Finds the positions where the maximum element occurs in the unified sequence.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2897 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, queso_require_not_equal_to_msg, UQ_FILE_EXTENSION_FOR_HDF_FORMAT, UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, UQ_FILE_EXTENSION_FOR_TXT_FORMAT, and QUESO::UQ_OK_RC.

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

2901 {
2902  queso_require_not_equal_to_msg(inputFileType, std::string(UQ_FILE_EXTENSION_FOR_TXT_FORMAT), std::string("reading txt files is not yet supported"));
2903  std::string fileType(inputFileType);
2904 #ifdef QUESO_HAS_HDF5
2905  // Do nothing
2906 #else
2907  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2908  if (m_env.subDisplayFile()) {
2909  *m_env.subDisplayFile() << "WARNING in ScalarSequence<T>::unifiedReadContents()"
2910  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2911  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2912  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2913  << "' instead..."
2914  << std::endl;
2915  }
2916  if (m_env.subRank() == 0) {
2917  std::cerr << "WARNING in ScalarSequence<T>::unifiedReadContents()"
2918  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2919  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2920  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2921  << "' instead..."
2922  << std::endl;
2923  }
2925  }
2926 #endif
2927 
2928  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2929  if (m_env.subDisplayFile()) {
2930  *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedReadContents()"
2931  << ": worldRank " << m_env.worldRank()
2932  << ", fullRank " << m_env.fullRank()
2933  << ", subEnvironment " << m_env.subId()
2934  << ", subRank " << m_env.subRank()
2935  << ", inter0Rank " << m_env.inter0Rank()
2936  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
2937  << ", fileName = " << fileName
2938  << ", subReadSize = " << subReadSize
2939  //<< ", unifiedReadSize = " << unifiedReadSize
2940  << std::endl;
2941  }
2942 
2943  this->resizeSequence(subReadSize);
2944 
2945  if (m_env.inter0Rank() >= 0) {
2946  double unifiedReadSize = subReadSize*m_env.inter0Comm().NumProc();
2947 
2948  // In the logic below, the id of a line' begins with value 0 (zero)
2949  unsigned int idOfMyFirstLine = 1 + m_env.inter0Rank()*subReadSize;
2950  unsigned int idOfMyLastLine = (1 + m_env.inter0Rank())*subReadSize;
2951  unsigned int numParams = 1; // this->vectorSizeLocal();
2952 
2953  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) { // "m or hdf"
2954  if (m_env.inter0Rank() == (int) r) {
2955  // My turn
2956  FilePtrSetStruct unifiedFilePtrSet;
2957  if (m_env.openUnifiedInputFile(fileName,
2958  fileType,
2959  unifiedFilePtrSet)) {
2960  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
2961  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
2962  if (r == 0) {
2963  // Read number of chain positions in the file by taking care of the first line,
2964  // which resembles something like 'variable_name = zeros(n_positions,m_params);'
2965  std::string tmpString;
2966 
2967  // Read 'variable name' string
2968  *unifiedFilePtrSet.ifsVar >> tmpString;
2969  //std::cout << "Just read '" << tmpString << "'" << std::endl;
2970 
2971  // Read '=' sign
2972  *unifiedFilePtrSet.ifsVar >> tmpString;
2973  //std::cout << "Just read '" << tmpString << "'" << std::endl;
2974  queso_require_equal_to_msg(tmpString, std::string("="), std::string("string should be the '=' sign"));
2975 
2976  // Read 'zeros(n_positions,n_params)' string
2977  *unifiedFilePtrSet.ifsVar >> tmpString;
2978  //std::cout << "Just read '" << tmpString << "'" << std::endl;
2979  unsigned int posInTmpString = 6;
2980 
2981  // Isolate 'n_positions' in a string
2982  //char nPositionsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
2983  std::string nPositionsString((size_t) (tmpString.size()-posInTmpString+1),' ');
2984  unsigned int posInPositionsString = 0;
2985  do {
2986  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ',' not found in first line of file");
2987  nPositionsString[posInPositionsString++] = tmpString[posInTmpString++];
2988  } while (tmpString[posInTmpString] != ',');
2989  nPositionsString[posInPositionsString] = '\0';
2990 
2991  // Isolate 'n_params' in a string
2992  posInTmpString++; // Avoid reading ',' char
2993  //char nParamsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
2994  std::string nParamsString((size_t) (tmpString.size()-posInTmpString+1),' ');
2995  unsigned int posInParamsString = 0;
2996  do {
2997  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ')' not found in first line of file");
2998  nParamsString[posInParamsString++] = tmpString[posInTmpString++];
2999  } while (tmpString[posInTmpString] != ')');
3000  nParamsString[posInParamsString] = '\0';
3001 
3002  // Convert 'n_positions' and 'n_params' strings to numbers
3003  unsigned int sizeOfChainInFile = (unsigned int) strtod(nPositionsString.c_str(),NULL);
3004  unsigned int numParamsInFile = (unsigned int) strtod(nParamsString.c_str(), NULL);
3005  if (m_env.subDisplayFile()) {
3006  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedReadContents()"
3007  << ": worldRank " << m_env.worldRank()
3008  << ", fullRank " << m_env.fullRank()
3009  << ", sizeOfChainInFile = " << sizeOfChainInFile
3010  << ", numParamsInFile = " << numParamsInFile
3011  << std::endl;
3012  }
3013 
3014  // Check if [size of chain in file] >= [requested unified sequence size]
3015  queso_require_greater_equal_msg(sizeOfChainInFile, unifiedReadSize, "size of chain in file is not big enough");
3016 
3017  // Check if [num params in file] == [num params in current chain]
3018  queso_require_equal_to_msg(numParamsInFile, numParams, "number of parameters of chain in file is different than number of parameters in this chain object");
3019  } // if (r == 0)
3020 
3021  // Code common to any core in 'inter0Comm', including core of rank 0
3022  unsigned int maxCharsPerLine = 64*numParams; // Up to about 60 characters to represent each parameter value
3023 
3024  unsigned int lineId = 0;
3025  while (lineId < idOfMyFirstLine) {
3026  unifiedFilePtrSet.ifsVar->ignore(maxCharsPerLine,'\n');
3027  lineId++;
3028  };
3029 
3030  if (r == 0) {
3031  // Take care of initial part of the first data line,
3032  // which resembles something like 'variable_name = [value1 value2 ...'
3033  std::string tmpString;
3034 
3035  // Read 'variable name' string
3036  *unifiedFilePtrSet.ifsVar >> tmpString;
3037  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
3038 
3039  // Read '=' sign
3040  *unifiedFilePtrSet.ifsVar >> tmpString;
3041  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
3042  queso_require_equal_to_msg(tmpString, std::string("="), std::string("in core 0, string should be the '=' sign"));
3043 
3044  // Take into account the ' [' portion
3045  std::streampos tmpPos = unifiedFilePtrSet.ifsVar->tellg();
3046  unifiedFilePtrSet.ifsVar->seekg(tmpPos+(std::streampos)2);
3047  }
3048 
3049  T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
3050  while (lineId <= idOfMyLastLine) {
3051  for (unsigned int i = 0; i < numParams; ++i) {
3052  *unifiedFilePtrSet.ifsVar >> tmpScalar;
3053  }
3054  m_seq[lineId - idOfMyFirstLine] = tmpScalar;
3055  lineId++;
3056  };
3057  }
3058 #ifdef QUESO_HAS_HDF5
3059  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
3060  if (r == 0) {
3061  hid_t dataset = H5Dopen2(unifiedFilePtrSet.h5Var,
3062  "data",
3063  H5P_DEFAULT); // Dataset access property list
3064  hid_t datatype = H5Dget_type(dataset);
3065  H5T_class_t t_class = H5Tget_class(datatype);
3066  queso_require_equal_to_msg(t_class, H5T_FLOAT, "t_class is not H5T_DOUBLE");
3067  hid_t dataspace = H5Dget_space(dataset);
3068  int rank = H5Sget_simple_extent_ndims(dataspace);
3069  queso_require_equal_to_msg(rank, 2, "hdf rank is not 2");
3070  hsize_t dims_in[2];
3071  int status_n;
3072  status_n = H5Sget_simple_extent_dims(dataspace, dims_in, NULL);
3073  if (status_n) {}; // just to remove compiler warning
3074  //std::cout << "In ScalarSequence<T>::unifiedReadContents()"
3075  // << ": dims_in[0] = " << dims_in[0]
3076  // << ", dims_in[1] = " << dims_in[1]
3077  // << std::endl;
3078  queso_require_equal_to_msg(dims_in[0], numParams, "dims_in[0] is not equal to 'numParams'");
3079  queso_require_greater_equal_msg(dims_in[1], subReadSize, "dims_in[1] is smaller that requested 'subReadSize'");
3080 
3081  struct timeval timevalBegin;
3082  int iRC = UQ_OK_RC;
3083  iRC = gettimeofday(&timevalBegin,NULL);
3084  if (iRC) {}; // just to remove compiler warning
3085 
3086  unsigned int chainSizeIn = (unsigned int) dims_in[1];
3087  //double* dataIn[numParams]; // avoid compiler warning
3088  std::vector<double*> dataIn((size_t) numParams,NULL);
3089  dataIn[0] = (double*) malloc(numParams*chainSizeIn*sizeof(double));
3090  for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1'
3091  dataIn[i] = dataIn[i-1] + chainSizeIn; // Yes, just 'chainSizeIn', not 'chainSizeIn*sizeof(double)'
3092  }
3093  //std::cout << "In ScalarSequence<T>::unifiedReadContents(): h5 case, memory allocated" << std::endl;
3094  herr_t status;
3095  status = H5Dread(dataset,
3096  H5T_NATIVE_DOUBLE,
3097  H5S_ALL,
3098  dataspace,
3099  H5P_DEFAULT,
3100  dataIn[0]);
3101  if (status) {}; // just to remove compiler warning
3102  //std::cout << "In ScalarSequence<T>::unifiedReadContents(): h5 case, data read" << std::endl;
3103  T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
3104  for (unsigned int j = 0; j < subReadSize; ++j) { // Yes, 'subReadSize', not 'chainSizeIn'
3105  for (unsigned int i = 0; i < numParams; ++i) {
3106  tmpScalar = dataIn[i][j];
3107  }
3108  m_seq[j] = tmpScalar;
3109  }
3110 
3111  double readTime = MiscGetEllapsedSeconds(&timevalBegin);
3112  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
3113  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedReadContents()"
3114  << ": worldRank " << m_env.worldRank()
3115  << ", fullRank " << m_env.fullRank()
3116  << ", subEnvironment " << m_env.subId()
3117  << ", subRank " << m_env.subRank()
3118  << ", inter0Rank " << m_env.inter0Rank()
3119  << ", fileName = " << fileName
3120  << ", numParams = " << numParams
3121  << ", chainSizeIn = " << chainSizeIn
3122  << ", subReadSize = " << subReadSize
3123  << ", readTime = " << readTime << " seconds"
3124  << std::endl;
3125  }
3126 
3127  H5Sclose(dataspace);
3128  H5Tclose(datatype);
3129  H5Dclose(dataset);
3130  //free(dataIn[0]); // related to the changes above for compiler warning
3131  for (unsigned int tmpIndex = 0; tmpIndex < dataIn.size(); tmpIndex++) {
3132  free (dataIn[tmpIndex]);
3133  }
3134  }
3135  else {
3136  queso_error_msg("hdf file type not supported for multiple sub-environments yet");
3137  }
3138  }
3139 #endif
3140  else {
3141  queso_error_msg("invalid file type");
3142  }
3143  m_env.closeFile(unifiedFilePtrSet,fileType);
3144  } // if (m_env.openUnifiedInputFile())
3145  } // if (m_env.inter0Rank() == (int) r)
3146  m_env.inter0Comm().Barrier();
3147  } // for r
3148  } // if (m_env.inter0Rank() >= 0)
3149  else {
3150  T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
3151  for (unsigned int i = 1; i < subReadSize; ++i) {
3152  m_seq[i] = tmpScalar;
3153  }
3154  }
3155 
3156  if (m_env.subDisplayFile()) {
3157  *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedReadContents()"
3158  << ", fileName = " << fileName
3159  << std::endl;
3160  }
3161  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
3162 
3163  return;
3164 }
int NumProc() const
Returns total number of processes.
Definition: MpiComm.C:123
const int UQ_OK_RC
Definition: Defines.h:89
int worldRank() const
Returns the process world rank.
Definition: Environment.C:262
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
Definition: MpiComm.C:164
double MiscGetEllapsedSeconds(struct timeval *timeval0)
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
int subRank() const
Access function for sub-rank.
Definition: Environment.C:287
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
#define UQ_FILE_EXTENSION_FOR_HDF_FORMAT
Definition: Defines.h:104
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:75
#define UQ_FILE_EXTENSION_FOR_TXT_FORMAT
Definition: Defines.h:103
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:341
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
Definition: Environment.C:1083
int fullRank() const
Returns the process full rank.
Definition: Environment.C:268
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
#define queso_require_not_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:74
istream * dataIn
Definition: ann_sample.cpp:58
unsigned int displayVerbosity() const
Definition: Environment.C:449
std::vector< T > m_seq
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:102
bool openUnifiedInputFile(const std::string &fileName, const std::string &fileType, FilePtrSetStruct &filePtrSet) const
Opens the unified input file.
Definition: Environment.C:990
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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 1158 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, and RawValue_MPI_SUM.

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

1163 {
1164  if (m_env.numSubEnvironments() == 1) {
1165  return this->subSampleStd(initialPos,
1166  numPos,
1167  unifiedMeanValue);
1168  }
1169 
1170  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1171 
1172  T unifiedStdValue = 0.;
1173  if (useOnlyInter0Comm) {
1174  if (m_env.inter0Rank() >= 0) {
1175  bool bRC = ((initialPos < this->subSequenceSize()) &&
1176  (0 < numPos ) &&
1177  ((initialPos+numPos) <= this->subSequenceSize()));
1178  queso_require_msg(bRC, "invalid input data");
1179 
1180  unsigned int finalPosPlus1 = initialPos + numPos;
1181  T diff;
1182  T localStdValue = 0.;
1183  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1184  diff = m_seq[j] - unifiedMeanValue;
1185  localStdValue += diff*diff;
1186  }
1187 
1188  unsigned int unifiedNumPos = 0;
1189  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
1190  "ScalarSequence<T>::unifiedSampleStd()",
1191  "failed MPI.Allreduce() for numPos");
1192 
1193  m_env.inter0Comm().template Allreduce<double>(&localStdValue, &unifiedStdValue, (int) 1, RawValue_MPI_SUM,
1194  "ScalarSequence<T>::unifiedSampleStd()",
1195  "failed MPI.Allreduce() for stdValue");
1196 
1197  unifiedStdValue /= (((T) unifiedNumPos) - 1.);
1198  unifiedStdValue = sqrt(unifiedStdValue);
1199  }
1200  else {
1201  // Node not in the 'inter0' communicator
1202  this->subSampleStd(initialPos,
1203  numPos,
1204  unifiedMeanValue);
1205  }
1206  }
1207  else {
1208  queso_error_msg("parallel vectors not supported yet");
1209  }
1210 
1211  //m_env.fullComm().Barrier();
1212 
1213  return unifiedStdValue;
1214 }
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
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...
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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
samVecis the vector of the calculated sample variance of the unified sequence of scalars.

Definition at line 1071 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, and RawValue_MPI_SUM.

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

1076 {
1077  if (m_env.numSubEnvironments() == 1) {
1078  return this->subSampleVarianceExtra(initialPos,
1079  numPos,
1080  unifiedMeanValue);
1081  }
1082 
1083  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1084 
1085  T unifiedSamValue = 0.;
1086  if (useOnlyInter0Comm) {
1087  if (m_env.inter0Rank() >= 0) {
1088  bool bRC = ((initialPos < this->subSequenceSize()) &&
1089  (0 < numPos ) &&
1090  ((initialPos+numPos) <= this->subSequenceSize()));
1091  queso_require_msg(bRC, "invalid input data");
1092 
1093  unsigned int finalPosPlus1 = initialPos + numPos;
1094  T diff;
1095  T localSamValue = 0.;
1096  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1097  diff = m_seq[j] - unifiedMeanValue;
1098  localSamValue += diff*diff;
1099  }
1100 
1101  unsigned int unifiedNumPos = 0;
1102  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
1103  "ScalarSequence<T>::unifiedSampleVarianceExtra()",
1104  "failed MPI.Allreduce() for numPos");
1105 
1106  m_env.inter0Comm().template Allreduce<double>(&localSamValue, &unifiedSamValue, (int) 1, RawValue_MPI_SUM,
1107  "ScalarSequence<T>::unifiedSampleVarianceExtra()",
1108  "failed MPI.Allreduce() for samValue");
1109 
1110  unifiedSamValue /= (((T) unifiedNumPos) - 1.);
1111  }
1112  else {
1113  // Node not in the 'inter0' communicator
1114  this->subSampleVarianceExtra(initialPos,
1115  numPos,
1116  unifiedMeanValue);
1117  }
1118  }
1119  else {
1120  queso_error_msg("parallel vectors not supported yet");
1121  }
1122 
1123  //m_env.fullComm().Barrier();
1124 
1125  return unifiedSamValue;
1126 }
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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 init...
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 455 of file ScalarSequence.C.

456 {
457  if (m_unifiedSampleVariancePlain == NULL) {
458  m_unifiedSampleVariancePlain = new T(0.);
460  }
461 
463 }
const T & unifiedMeanPlain(bool useOnlyInter0Comm) const
Finds the mean value of the unified sequence of scalars.
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 ...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2239 of file ScalarSequence.C.

References queso_error_msg, queso_require_msg, and RawValue_MPI_SUM.

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

2244 {
2245  if (m_env.numSubEnvironments() == 1) {
2246  return this->subScaleForKde(initialPos,
2247  unifiedIqrValue,
2248  kdeDimension);
2249  }
2250 
2251  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
2252 
2253  T unifiedScaleValue = 0.;
2254  if (useOnlyInter0Comm) {
2255  if (m_env.inter0Rank() >= 0) {
2256  bool bRC = (initialPos < this->subSequenceSize());
2257  queso_require_msg(bRC, "invalid input data");
2258 
2259  unsigned int localDataSize = this->subSequenceSize() - initialPos;
2260 
2261  T unifiedMeanValue = this->unifiedMeanExtra(useOnlyInter0Comm,
2262  initialPos,
2263  localDataSize);
2264 
2265  T unifiedSamValue = this->unifiedSampleVarianceExtra(useOnlyInter0Comm,
2266  initialPos,
2267  localDataSize,
2268  unifiedMeanValue);
2269 
2270  unsigned int unifiedDataSize = 0;
2271  m_env.inter0Comm().template Allreduce<unsigned int>(&localDataSize, &unifiedDataSize, (int) 1, RawValue_MPI_SUM,
2272  "ScalarSequence<T>::unifiedScaleForKde()",
2273  "failed MPI.Allreduce() for data size");
2274 
2275  if (unifiedIqrValue <= 0.) {
2276  unifiedScaleValue = 1.06*std::sqrt(unifiedSamValue)/std::pow(unifiedDataSize,1./(4. + ((double) kdeDimension)));
2277  }
2278  else {
2279  unifiedScaleValue = 1.06*std::min(std::sqrt(unifiedSamValue),unifiedIqrValue/1.34)/std::pow(unifiedDataSize,1./(4. + ((double) kdeDimension)));
2280  }
2281 
2282  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2283  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedScaleForKde()"
2284  << ": unifiedIqrValue = " << unifiedIqrValue
2285  << ", unifiedMeanValue = " << unifiedMeanValue
2286  << ", unifiedSamValue = " << unifiedSamValue
2287  << ", unifiedDataSize = " << unifiedDataSize
2288  << ", unifiedScaleValue = " << unifiedScaleValue
2289  << std::endl;
2290  }
2291  }
2292  else {
2293  // Node not in the 'inter0' communicator
2294  unifiedScaleValue = this->subScaleForKde(initialPos,
2295  unifiedIqrValue,
2296  kdeDimension);
2297  }
2298  }
2299  else {
2300  queso_error_msg("parallel vectors not supported yet");
2301  }
2302 
2303  //m_env.fullComm().Barrier();
2304 
2305  return unifiedScaleValue;
2306 }
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
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-sequenc...
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
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 ...
unsigned int displayVerbosity() const
Definition: Environment.C:449
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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...
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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, and RawValue_MPI_SUM.

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().

144 {
145  if (m_env.numSubEnvironments() == 1) {
146  return this->subSequenceSize();
147  }
148 
149  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
150 
151  unsigned int unifiedNumSamples = 0;
152  if (useOnlyInter0Comm) {
153  if (m_env.inter0Rank() >= 0) {
154  unsigned int subNumSamples = this->subSequenceSize();
155  m_env.inter0Comm().template Allreduce<unsigned int>(&subNumSamples, &unifiedNumSamples, (int) 1, RawValue_MPI_SUM,
156  "ScalarSequence<T>::unifiedSequenceSize()",
157  "failed MPI.Allreduce() for unifiedSequenceSize()");
158  }
159  else {
160  // Node not in the 'inter0' communicator
161  unifiedNumSamples = this->subSequenceSize();
162  }
163  }
164  else {
165  queso_error_msg("parallel vectors not supported yet");
166  }
167 
168  return unifiedNumSamples;
169 }
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
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 1891 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().

1895 {
1896  if (m_env.numSubEnvironments() == 1) {
1897  return this->subSort(initialPos,unifiedSortedSequence);
1898  }
1899 
1900  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
1901 
1902  if (useOnlyInter0Comm) {
1903  if (m_env.inter0Rank() >= 0) {
1904  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), beginning logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
1905 
1906  unsigned int localNumPos = this->subSequenceSize() - initialPos;
1907 
1908  std::vector<T> leafData(localNumPos,0.);
1909  this->extractRawData(0,
1910  1,
1911  localNumPos,
1912  leafData);
1913 
1914  if (m_env.inter0Rank() == 0) {
1915  int minus1NumTreeLevels = 0;
1916  int power2NumTreeNodes = 1;
1917 
1918  while (power2NumTreeNodes < m_env.inter0Comm().NumProc()) {
1919  power2NumTreeNodes += power2NumTreeNodes;
1920  minus1NumTreeLevels++;
1921  }
1922 
1923  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1924  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedSort()"
1925  << ": sorting tree has " << m_env.inter0Comm().NumProc()
1926  << " nodes and " << minus1NumTreeLevels+1
1927  << " levels"
1928  << std::endl;
1929  }
1930 
1931  this->parallelMerge(unifiedSortedSequence.rawData(),
1932  leafData,
1933  minus1NumTreeLevels);
1934  }
1935  else if (m_env.inter0Rank() > 0) { // KAUST
1936  unsigned int uintBuffer[1];
1937  RawType_MPI_Status status;
1939  "ScalarSequence<T>::unifiedSort()",
1940  "failed MPI.Recv() for init");
1941  //if (status) {}; // just to remove compiler warning
1942 
1943  unsigned int treeLevel = uintBuffer[0];
1944  this->parallelMerge(unifiedSortedSequence.rawData(),
1945  leafData,
1946  treeLevel);
1947  }
1948 
1949  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), returned from parallelMerge()",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
1950 
1951  // Broadcast
1952  unsigned int unifiedDataSize = unifiedSortedSequence.subSequenceSize();
1953  m_env.inter0Comm().Bcast((void *) &unifiedDataSize, (int) 1, RawValue_MPI_UNSIGNED, 0,
1954  "ScalarSequence<T>::unifiedSort()",
1955  "failed MPI.Bcast() for unified data size");
1956 
1957  unsigned int sumOfNumPos = 0;
1958  m_env.inter0Comm().template Allreduce<unsigned int>(&localNumPos, &sumOfNumPos, (int) 1, RawValue_MPI_SUM,
1959  "ScalarSequence<T>::unifiedSort()",
1960  "failed MPI.Allreduce() for data size");
1961 
1962  queso_require_equal_to_msg(sumOfNumPos, unifiedDataSize, "incompatible unified sizes");
1963 
1964  unifiedSortedSequence.resizeSequence(unifiedDataSize);
1965  m_env.inter0Comm().Bcast((void *) &unifiedSortedSequence.rawData()[0], (int) unifiedDataSize, RawValue_MPI_DOUBLE, 0,
1966  "ScalarSequence<T>::unifiedSort()",
1967  "failed MPI.Bcast() for unified data");
1968 
1969  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1970  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
1971  << ": tree node " << m_env.inter0Rank()
1972  << ", unifiedSortedSequence[0] = " << unifiedSortedSequence[0]
1973  << ", unifiedSortedSequence[" << unifiedSortedSequence.subSequenceSize()-1 << "] = " << unifiedSortedSequence[unifiedSortedSequence.subSequenceSize()-1]
1974  << std::endl;
1975  }
1976 
1977  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), ending logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
1978  }
1979  else {
1980  // Node not in the 'inter0' communicator
1981  this->subSort(initialPos,unifiedSortedSequence);
1982  }
1983  }
1984  else {
1985  queso_error_msg("parallel vectors not supported yet");
1986  }
1987 
1988  //m_env.fullComm().Barrier();
1989 
1990  return;
1991 }
int NumProc() const
Returns total number of processes.
Definition: MpiComm.C:123
int RawType_MPI_Status
Definition: MpiComm.h:62
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, std::vector< double > &rawData) const
Extracts the raw data.
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
void parallelMerge(std::vector< T > &sortedBuffer, const std::vector< T > &leafData, unsigned int treeLevel) const
Sorts/merges data in parallel using MPI.
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:67
#define RawValue_MPI_UNSIGNED
Definition: MpiComm.h:68
#define RawValue_MPI_ANY_SOURCE
Definition: MpiComm.h:64
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
#define SCALAR_SEQUENCE_INIT_MPI_MSG
unsigned int displayVerbosity() const
Definition: Environment.C:449
void Recv(void *buf, int count, RawType_MPI_Datatype datatype, int source, int tag, RawType_MPI_Status *status, const char *whereMsg, const char *whatMsg) const
Blocking receive of data from this process to another process.
Definition: MpiComm.C:297
void subSort()
Sorts the sequence of scalars in the private attribute m_seq.
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
void Bcast(void *buffer, int count, RawType_MPI_Datatype datatype, int root, const char *whereMsg, const char *whatMsg) const
Broadcast values from the root process to the slave processes.
Definition: MpiComm.C:181
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 613 of file ScalarSequence.C.

References queso_error_msg.

619 {
620  if (m_env.numSubEnvironments() == 1) {
621  return this->subUniformlySampledCdf(numEvaluationPoints,
622  unifiedMinDomainValue,
623  unifiedMaxDomainValue,
624  unifiedCdfValues);
625  }
626 
627  // KAUST2
628  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
629 
630  if (useOnlyInter0Comm) {
631  if (m_env.inter0Rank() >= 0) {
632  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
633  *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedUniformlySampledCdf()"
634  << std::endl;
635  }
636 
637  T unifiedTmpMinValue;
638  T unifiedTmpMaxValue;
639  std::vector<T> unifiedCenters(numEvaluationPoints,0.);
640  std::vector<unsigned int> unifiedBins (numEvaluationPoints,0);
641 
642  this->unifiedMinMaxExtra(useOnlyInter0Comm,
643  0, // initialPos
644  this->subSequenceSize(),
645  unifiedTmpMinValue,
646  unifiedTmpMaxValue);
647  this->unifiedHistogram(useOnlyInter0Comm,
648  0, // initialPos
649  unifiedTmpMinValue,
650  unifiedTmpMaxValue,
651  unifiedCenters,
652  unifiedBins);
653 
654  unifiedMinDomainValue = unifiedCenters[0];
655  unifiedMaxDomainValue = unifiedCenters[unifiedCenters.size()-1];
656 
657  unsigned int unifiedTotalSumOfBins = 0;
658  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
659  unifiedTotalSumOfBins += unifiedBins[i];
660  }
661 
662  std::vector<unsigned int> unifiedPartialSumsOfBins(numEvaluationPoints,0);
663  unifiedPartialSumsOfBins[0] = unifiedBins[0];
664  for (unsigned int i = 1; i < numEvaluationPoints; ++i) { // Yes, from '1'
665  unifiedPartialSumsOfBins[i] = unifiedPartialSumsOfBins[i-1] + unifiedBins[i];
666  }
667 
668  unifiedCdfValues.clear();
669  unifiedCdfValues.resize(numEvaluationPoints);
670  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
671  unifiedCdfValues[i] = ((T) unifiedPartialSumsOfBins[i])/((T) unifiedTotalSumOfBins);
672  }
673 
674  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
675  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
676  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedUniformlySampledCdf()"
677  << ": i = " << i
678  << ", unifiedTmpMinValue = " << unifiedTmpMinValue
679  << ", unifiedTmpMaxValue = " << unifiedTmpMaxValue
680  << ", unifiedBins = " << unifiedBins[i]
681  << ", unifiedCdfValue = " << unifiedCdfValues[i]
682  << ", unifiedPartialSumsOfBins = " << unifiedPartialSumsOfBins[i]
683  << ", unifiedTotalSumOfBins = " << unifiedTotalSumOfBins
684  << std::endl;
685  }
686  }
687 
688  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
689  *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedUniformlySampledCdf()"
690  << std::endl;
691  }
692  }
693  else {
694  // Node not in the 'inter0' communicator
695  this->subUniformlySampledCdf(numEvaluationPoints,
696  unifiedMinDomainValue,
697  unifiedMaxDomainValue,
698  unifiedCdfValues);
699  }
700  }
701  else {
702  queso_error_msg("parallel vectors not supported yet");
703  }
704 
705  //m_env.fullComm().Barrier();
706 
707  return;
708 }
void subUniformlySampledCdf(unsigned int numIntervals, T &minDomainValue, T &maxDomainValue, std::vector< T > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
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 starti...
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
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.
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
unsigned int displayVerbosity() const
Definition: Environment.C:449
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2672 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, UQ_FILE_EXTENSION_FOR_TXT_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().

2675 {
2676  std::string fileType(inputFileType);
2677 #ifdef QUESO_HAS_HDF5
2678  // Do nothing
2679 #else
2680  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2681  if (m_env.subDisplayFile()) {
2682  *m_env.subDisplayFile() << "WARNING in ScalarSequence<T>::unifiedWriteContents()"
2683  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2684  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2685  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2686  << "' instead..."
2687  << std::endl;
2688  }
2689  if (m_env.subRank() == 0) {
2690  std::cerr << "WARNING in ScalarSequence<T>::unifiedWriteContents()"
2691  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2692  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2693  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2694  << "' instead..."
2695  << std::endl;
2696  }
2698  }
2699 #endif
2700 
2701  // All processors in 'fullComm' should call this routine...
2702 
2703  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2704  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
2705  *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedWriteContents()"
2706  << ": worldRank " << m_env.worldRank()
2707  << ", subEnvironment " << m_env.subId()
2708  << ", subRank " << m_env.subRank()
2709  << ", inter0Rank " << m_env.inter0Rank()
2710  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
2711  << ", fileName = " << fileName
2712  << ", fileType = " << fileType
2713  << std::endl;
2714  }
2715 
2716  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
2717 
2718  if (m_env.inter0Rank() >= 0) {
2719  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) {
2720  if (m_env.inter0Rank() == (int) r) {
2721  // My turn
2722  FilePtrSetStruct unifiedFilePtrSet;
2723  // bool writeOver = (r == 0);
2724  bool writeOver = false; // A 'true' causes problems when the user chooses (via options
2725  // in the input file) to use just one file for all outputs.
2726  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 000 \n" << std::endl;
2727  if (m_env.openUnifiedOutputFile(fileName,
2728  fileType, // "m or hdf"
2729  writeOver,
2730  unifiedFilePtrSet)) {
2731  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 001 \n" << std::endl;
2732 
2733  unsigned int chainSize = this->subSequenceSize();
2734  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
2735  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
2736 
2737  // Write header info
2738  if (r == 0) {
2739  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2740  this->writeUnifiedMatlabHeader(*unifiedFilePtrSet.ofsVar,
2741  this->subSequenceSize()*m_env.inter0Comm().NumProc());
2742  }
2743  else { // It's definitely txt if we get in here
2744  this->writeTxtHeader(*unifiedFilePtrSet.ofsVar,
2745  this->subSequenceSize()*m_env.inter0Comm().NumProc());
2746  }
2747  }
2748 
2749  for (unsigned int j = 0; j < chainSize; ++j) {
2750  *unifiedFilePtrSet.ofsVar << m_seq[j]
2751  << std::endl;
2752  }
2753 
2754  m_env.closeFile(unifiedFilePtrSet,fileType);
2755  }
2756 #ifdef QUESO_HAS_HDF5
2757  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2758  unsigned int numParams = 1; // m_vectorSpace.dimLocal();
2759  if (r == 0) {
2760  hid_t datatype = H5Tcopy(H5T_NATIVE_DOUBLE);
2761  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data type created" << std::endl;
2762  hsize_t dimsf[1];
2763  dimsf[0] = chainSize;
2764  hid_t dataspace = H5Screate_simple(1, dimsf, NULL); // HDF5_rank = 2
2765  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data space created" << std::endl;
2766  hid_t dataset = H5Dcreate2(unifiedFilePtrSet.h5Var,
2767  "data",
2768  datatype,
2769  dataspace,
2770  H5P_DEFAULT, // Link creation property list
2771  H5P_DEFAULT, // Dataset creation property list
2772  H5P_DEFAULT); // Dataset access property list
2773  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data set created" << std::endl;
2774 
2775  struct timeval timevalBegin;
2776  int iRC = UQ_OK_RC;
2777  iRC = gettimeofday(&timevalBegin,NULL);
2778  if (iRC) {}; // just to remove compiler warning
2779 
2780  herr_t status;
2781  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 002 \n" << std::endl;
2782  status = H5Dwrite(dataset,
2783  H5T_NATIVE_DOUBLE,
2784  H5S_ALL,
2785  H5S_ALL,
2786  H5P_DEFAULT,
2787  &m_seq[0]);
2788  if (status) {}; // just to remove compiler warning
2789 
2790  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 003 \n" << std::endl;
2791  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data written" << std::endl;
2792 
2793  double writeTime = MiscGetEllapsedSeconds(&timevalBegin);
2794  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2795  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedWriteContents()"
2796  << ": worldRank " << m_env.worldRank()
2797  << ", fullRank " << m_env.fullRank()
2798  << ", subEnvironment " << m_env.subId()
2799  << ", subRank " << m_env.subRank()
2800  << ", inter0Rank " << m_env.inter0Rank()
2801  << ", fileName = " << fileName
2802  << ", numParams = " << numParams
2803  << ", chainSize = " << chainSize
2804  << ", writeTime = " << writeTime << " seconds"
2805  << std::endl;
2806  }
2807 
2808  H5Dclose(dataset);
2809  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data set closed" << std::endl;
2810  H5Sclose(dataspace);
2811  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data space closed" << std::endl;
2812  H5Tclose(datatype);
2813  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data type closed" << std::endl;
2814  }
2815  else {
2816  queso_error_msg("hdf file type not supported for multiple sub-environments yet");
2817  }
2818  }
2819 #endif
2820  } // if (m_env.openUnifiedOutputFile())
2821  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 004 \n" << std::endl;
2822  } // if (m_env.inter0Rank() == (int) r)
2823  m_env.inter0Comm().Barrier();
2824  } // for r
2825 
2826  if (m_env.inter0Rank() == 0) {
2827  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
2828  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
2829  FilePtrSetStruct unifiedFilePtrSet;
2830  if (m_env.openUnifiedOutputFile(fileName,
2831  fileType,
2832  false, // Yes, 'writeOver = false' in order to close the array for matlab
2833  unifiedFilePtrSet)) {
2834 
2835  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2836  *unifiedFilePtrSet.ofsVar << "];\n";
2837  }
2838 
2839  m_env.closeFile(unifiedFilePtrSet,fileType);
2840  }
2841  }
2842  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2843  // Do nothing
2844  }
2845  else {
2846  queso_error_msg("invalid file type");
2847  }
2848  }
2849  } // if (m_env.inter0Rank() >= 0)
2850 
2851  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
2852  *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedWriteContents()"
2853  << ", fileName = " << fileName
2854  << ", fileType = " << fileType
2855  << std::endl;
2856  }
2857  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2858 
2859  return;
2860 }
int NumProc() const
Returns total number of processes.
Definition: MpiComm.C:123
const int UQ_OK_RC
Definition: Defines.h:89
int worldRank() const
Returns the process world rank.
Definition: Environment.C:262
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
Definition: MpiComm.C:164
double MiscGetEllapsedSeconds(struct timeval *timeval0)
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
int subRank() const
Access function for sub-rank.
Definition: Environment.C:287
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
#define UQ_FILE_EXTENSION_FOR_HDF_FORMAT
Definition: Defines.h:104
#define UQ_FILE_EXTENSION_FOR_TXT_FORMAT
Definition: Defines.h:103
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:341
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
Definition: Environment.C:1083
int fullRank() const
Returns the process full rank.
Definition: Environment.C:268
unsigned int displayVerbosity() const
Definition: Environment.C:449
std::vector< T > m_seq
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:102
void writeUnifiedMatlabHeader(std::ofstream &ofs, double sequenceSize) const
Helper function to write header info for matlab files from all chains.
#define queso_error_msg(msg)
Definition: asserts.h:47
void writeTxtHeader(std::ofstream &ofs, double sequenceSize) const
Helper function to write txt info for matlab files.
const BaseEnvironment & m_env
bool openUnifiedOutputFile(const std::string &fileName, const std::string &fileType, bool writeOver, FilePtrSetStruct &filePtrSet) const
Opens a unified output file, that will contain data from all sub-environments.
Definition: Environment.C:725
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
template<class T >
void QUESO::ScalarSequence< T >::writeSubMatlabHeader ( std::ofstream &  ofs,
double  sequenceSize 
) const
private

Helper function to write header info for matlab files from one chain.

Definition at line 2876 of file ScalarSequence.C.

2878 {
2879  ofs << m_name << "_sub" << m_env.subIdString() << " = zeros(" << sequenceSize
2880  << "," << 1
2881  << ");"
2882  << std::endl;
2883  ofs << m_name << "_sub" << m_env.subIdString() << " = [";
2884 }
const std::string & subIdString() const
Access to the attribute m_subIdString; which stores the string for the sub-environment, and it will be used, for instance, to create the output files for each sub-environment.
Definition: Environment.C:347
const BaseEnvironment & m_env
template<class T >
void QUESO::ScalarSequence< T >::writeTxtHeader ( std::ofstream &  ofs,
double  sequenceSize 
) const
private

Helper function to write txt info for matlab files.

Definition at line 2888 of file ScalarSequence.C.

2890 {
2891  ofs << sequenceSize << " " << 1 << std::endl;
2892 }
template<class T >
void QUESO::ScalarSequence< T >::writeUnifiedMatlabHeader ( std::ofstream &  ofs,
double  sequenceSize 
) const
private

Helper function to write header info for matlab files from all chains.

Definition at line 2864 of file ScalarSequence.C.

2866 {
2867  ofs << m_name << "_unified" << " = zeros(" << sequenceSize
2868  << "," << 1
2869  << ");"
2870  << std::endl;
2871  ofs << m_name << "_unified" << " = [";
2872 }

Member Data Documentation

template<class T = double>
const BaseEnvironment& QUESO::ScalarSequence< T >::m_env
private

Definition at line 525 of file ScalarSequence.h.

template<class T = double>
std::string QUESO::ScalarSequence< T >::m_name
private

Definition at line 526 of file ScalarSequence.h.

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

template<class T = double>
std::vector<T> QUESO::ScalarSequence< T >::m_seq
private
template<class T = double>
T* QUESO::ScalarSequence< T >::m_subMaxPlain
mutableprivate

Definition at line 531 of file ScalarSequence.h.

template<class T = double>
T* QUESO::ScalarSequence< T >::m_subMeanPlain
mutableprivate

Definition at line 533 of file ScalarSequence.h.

template<class T = double>
T* QUESO::ScalarSequence< T >::m_subMedianPlain
mutableprivate

Definition at line 535 of file ScalarSequence.h.

template<class T = double>
T* QUESO::ScalarSequence< T >::m_subMinPlain
mutableprivate

Definition at line 529 of file ScalarSequence.h.

template<class T = double>
T* QUESO::ScalarSequence< T >::m_subSampleVariancePlain
mutableprivate

Definition at line 537 of file ScalarSequence.h.

template<class T = double>
T* QUESO::ScalarSequence< T >::m_unifiedMaxPlain
mutableprivate

Definition at line 532 of file ScalarSequence.h.

template<class T = double>
T* QUESO::ScalarSequence< T >::m_unifiedMeanPlain
mutableprivate

Definition at line 534 of file ScalarSequence.h.

template<class T = double>
T* QUESO::ScalarSequence< T >::m_unifiedMedianPlain
mutableprivate

Definition at line 536 of file ScalarSequence.h.

template<class T = double>
T* QUESO::ScalarSequence< T >::m_unifiedMinPlain
mutableprivate

Definition at line 530 of file ScalarSequence.h.

template<class T = double>
T* QUESO::ScalarSequence< T >::m_unifiedSampleVariancePlain
mutableprivate

Definition at line 538 of file ScalarSequence.h.


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

Generated on Tue Nov 29 2016 10:53:13 for queso-0.56.0 by  doxygen 1.8.5