queso-0.55.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 & m_env
const BaseEnvironment & env() const
Access to QUESO environment.
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 2499 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().

2503 {
2504  queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+1), "srcInitialPos is too big");
2505 
2506  queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+srcNumPos), "srcNumPos is too big");
2507 
2509  unsigned int currentSize = this->subSequenceSize();
2510  m_seq.resize(currentSize+srcNumPos,0.);
2511  for (unsigned int i = 0; i < srcNumPos; ++i) {
2512  m_seq[currentSize+i] = src.m_seq[srcInitialPos+i];
2513  }
2514 
2515  return;
2516 }
std::vector< T > m_seq
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:90
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
void deleteStoredScalars()
Deletes all stored 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 }
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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.
T autoCovariance(unsigned int initialPos, unsigned int numPos, const T &meanValue, unsigned int lag) const
Calculates the autocovariance.
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  unsigned int fftSize = 0;
1370  {
1371 #warning WTF are 4 lines of unused code doing here? - RHS
1372  double tmp = log((double) maxLag)/log(2.);
1373  double fractionalPart = tmp - ((double) ((unsigned int) tmp));
1374  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
1375  unsigned int fftSize1 = (unsigned int) std::pow(2.,tmp+1.); // Yes, tmp+1
1376 
1377  tmp = log((double) numPos)/log(2.);
1378  fractionalPart = tmp - ((double) ((unsigned int) tmp));
1379  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
1380  unsigned int fftSize2 = (unsigned int) std::pow(2.,tmp+1);
1381 
1382  fftSize = fftSize2;
1383  }
1384 
1385  std::vector<double> rawDataVec(numPos,0.);
1386  std::vector<std::complex<double> > resultData(0,std::complex<double>(0.,0.));
1387  Fft<T> fftObj(m_env);
1388 
1389  // Forward FFT
1390  this->extractRawData(initialPos,
1391  1, // spacing
1392  numPos,
1393  rawDataVec);
1394  T meanValue = this->subMeanExtra(initialPos,
1395  numPos);
1396  for (unsigned int j = 0; j < numPos; ++j) {
1397  rawDataVec[j] -= meanValue; // IMPORTANT
1398  }
1399 
1400  rawDataVec.resize(fftSize,0.);
1401 
1402  //if (m_env.subDisplayFile()) {
1403  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1404  // << ": about to call fftObj.forward()"
1405  // << " with rawDataVec.size() = " << rawDataVec.size()
1406  // << ", fftSize = " << fftSize
1407  // << ", resultData.size() = " << resultData.size()
1408  // << std::endl;
1409  //}
1410  fftObj.forward(rawDataVec,fftSize,resultData);
1411 
1412  // Inverse FFT
1413  for (unsigned int j = 0; j < fftSize; ++j) {
1414  rawDataVec[j] = std::norm(resultData[j]);
1415  }
1416  //if (m_env.subDisplayFile()) {
1417  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1418  // << ": about to call fftObj.inverse()"
1419  // << " with rawDataVec.size() = " << rawDataVec.size()
1420  // << ", fftSize = " << fftSize
1421  // << ", resultData.size() = " << resultData.size()
1422  // << std::endl;
1423  //}
1424  fftObj.inverse(rawDataVec,fftSize,resultData);
1425  //if (m_env.subDisplayFile()) {
1426  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1427  // << ": returned succesfully from fftObj.inverse()"
1428  // << std::endl;
1429  //}
1430 
1431  // Prepare return data
1432  autoCorrs.resize(maxLag+1,0.); // Yes, +1
1433  for (unsigned int j = 0; j < autoCorrs.size(); ++j) {
1434  double ratio = ((double) j)/((double) (numPos-1));
1435  autoCorrs[j] = ( resultData[j].real()/resultData[0].real() )*(1.-ratio);
1436  }
1437 
1438  return;
1439 }
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, std::vector< double > &rawData) const
Extracts the raw data.
const BaseEnvironment & m_env
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
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 1443 of file ScalarSequence.C.

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

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

References queso_error_msg, and queso_not_implemented.

2474 {
2475  double resultValue = 0.;
2476 
2477  // FIX ME: put logic if (numSubEnvs == 1) ...
2478 
2479  if (useOnlyInter0Comm) {
2480  if (m_env.inter0Rank() >= 0) {
2482  }
2483  else {
2484  // Node not in the 'inter0' communicator
2485  // Do nothing
2486  }
2487  }
2488  else {
2489  queso_error_msg("parallel vectors not supported yet");
2490  }
2491 
2492  //m_env.fullComm().Barrier();
2493 
2494  return resultValue;
2495 }
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:265
#define queso_not_implemented()
Definition: asserts.h:56
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 resetValues(unsigned int initialPos, unsigned int)
Sets numPos values of the sequence to zero, starting at position initialPos.
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
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 3179 of file ScalarSequence.C.

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

3180 {
3181  m_name = src.m_name;
3182  m_seq.clear();
3183  m_seq.resize(src.subSequenceSize(),0.);
3184  for (unsigned int i = 0; i < m_seq.size(); ++i) {
3185  m_seq[i] = src.m_seq[i];
3186  }
3188 
3189  return;
3190 }
std::vector< T > m_seq
void deleteStoredScalars()
Deletes all stored scalars.
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 }
std::vector< T > m_seq
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
std::vector< T >::iterator seqScalarPositionIteratorTypedef
void deleteStoredScalars()
Deletes all stored 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 3235 of file ScalarSequence.C.

3240 {
3241  rawDataVec.resize(numPos);
3242  if (spacing == 1) {
3243  for (unsigned int j = 0; j < numPos; ++j) {
3244  rawDataVec[j] = m_seq[initialPos+j ];
3245  }
3246  }
3247  else {
3248  for (unsigned int j = 0; j < numPos; ++j) {
3249  rawDataVec[j] = m_seq[initialPos+j*spacing];
3250  }
3251  }
3252 
3253  return;
3254 }
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 3194 of file ScalarSequence.C.

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

3199 {
3200  scalarSeq.resizeSequence(numPos);
3201  if (spacing == 1) {
3202  for (unsigned int j = 0; j < numPos; ++j) {
3203  //if ((initialPos+j*spacing) > m_seq.size()) {
3204  // std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
3205  // << ": initialPos = " << initialPos
3206  // << ", spacing = " << spacing
3207  // << ", numPos = " << numPos
3208  // << ", j = " << j
3209  // << ", position got too large"
3210  // << std::endl;
3211  //}
3212  scalarSeq[j] = m_seq[initialPos+j ];
3213  }
3214  }
3215  else {
3216  for (unsigned int j = 0; j < numPos; ++j) {
3217  //if ((initialPos+j*spacing) > m_seq.size()) {
3218  // std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
3219  // << ": initialPos = " << initialPos
3220  // << ", spacing = " << spacing
3221  // << ", numPos = " << numPos
3222  // << ", j = " << j
3223  // << ", position got too large"
3224  // << std::endl;
3225  //}
3226  scalarSeq[j] = m_seq[initialPos+j*spacing];
3227  }
3228  }
3229 
3230  return;
3231 }
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 2431 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().

2434 {
2435  if (m_env.subDisplayFile()) {
2436  *m_env.subDisplayFile() << "Entering ScalarSequence<V,M>::filter()"
2437  << ": initialPos = " << initialPos
2438  << ", spacing = " << spacing
2439  << ", subSequenceSize = " << this->subSequenceSize()
2440  << std::endl;
2441  }
2442 
2443  unsigned int i = 0;
2444  unsigned int j = initialPos;
2445  unsigned int originalSubSequenceSize = this->subSequenceSize();
2446  while (j < originalSubSequenceSize) {
2447  if (i != j) {
2448  //*m_env.subDisplayFile() << i << "--" << j << " ";
2449  m_seq[i] = m_seq[j];
2450  }
2451  i++;
2452  j += spacing;
2453  }
2454 
2455  this->resizeSequence(i);
2456 
2457  if (m_env.subDisplayFile()) {
2458  *m_env.subDisplayFile() << "Leaving ScalarSequence<V,M>::filter()"
2459  << ": initialPos = " << initialPos
2460  << ", spacing = " << spacing
2461  << ", subSequenceSize = " << this->subSequenceSize()
2462  << std::endl;
2463  }
2464 
2465  return;
2466 }
std::vector< T > m_seq
const BaseEnvironment & m_env
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
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 }
unsigned int displayVerbosity() const
Definition: Environment.C:400
#define queso_error_msg(msg)
Definition: asserts.h:47
std::vector< T > m_seq
const BaseEnvironment & m_env
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:265
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
int NumProc() const
Returns total number of processes.
Definition: MpiComm.C:123
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:271
unsigned int subSequenceSize() const
Size of the sub-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
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:67
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:87
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:87
std::vector< T > m_seq
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
void deleteStoredScalars()
Deletes all stored 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 3276 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.

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

The sequence of scalars. Access to private attribute m_seq.

Definition at line 3258 of file ScalarSequence.C.

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

3259 {
3260  return m_seq;
3261 }
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 }
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.
void deleteStoredScalars()
Deletes all stored 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 }
std::vector< T > m_seq
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
void deleteStoredScalars()
Deletes all stored 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 }
virtual double gaussianSample(double stdDev) const =0
Samples a value from a Gaussian distribution with standard deviation given by stdDev.
std::vector< T > m_seq
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
const RngBase * rngObject() const
Access to the RNG object.
Definition: Environment.C:421
void deleteStoredScalars()
Deletes all stored 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 }
std::vector< T > m_seq
const BaseEnvironment & m_env
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
virtual double uniformSample() const =0
Samples a value from a uniform distribution.
const RngBase * rngObject() const
Access to the RNG object.
Definition: Environment.C:421
void deleteStoredScalars()
Deletes all stored 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.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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...
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 1749 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

1755 {
1756  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1757 
1758  for (unsigned int j = 0; j < bins.size(); ++j) {
1759  bins[j] = 0;
1760  }
1761 
1762  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1763  double minCenter = minHorizontalValue - horizontalDelta/2.;
1764  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1765  gridValues = new UniformOneDGrid<T>(m_env,
1766  "",
1767  bins.size(),
1768  minCenter,
1769  maxCenter);
1770 
1771  unsigned int dataSize = this->subSequenceSize();
1772  for (unsigned int j = 0; j < dataSize; ++j) {
1773  double value = m_seq[j];
1774  if (value < minHorizontalValue) {
1775  bins[0] += value;
1776  }
1777  else if (value >= maxHorizontalValue) {
1778  bins[bins.size()-1] += value;
1779  }
1780  else {
1781  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1782  bins[index] += value;
1783  }
1784  }
1785 
1786  return;
1787 }
std::vector< T > m_seq
const BaseEnvironment & m_env
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:90
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 2321 of file ScalarSequence.C.

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

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

2326 {
2327  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2328  (0 < evaluationPositions.size()) &&
2329  (evaluationPositions.size() == densityValues.size() ));
2330  queso_require_msg(bRC, "invalid input data");
2331 
2332  unsigned int dataSize = this->subSequenceSize() - initialPos;
2333  unsigned int numEvals = evaluationPositions.size();
2334 
2335  double scaleInv = 1./scaleValue;
2336  for (unsigned int j = 0; j < numEvals; ++j) {
2337  double x = evaluationPositions[j];
2338  double value = 0.;
2339  for (unsigned int k = 0; k < dataSize; ++k) {
2340  double xk = m_seq[initialPos+k];
2341  value += MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
2342  }
2343  densityValues[j] = scaleInv * (value/(double) dataSize);
2344  }
2345 
2346  return;
2347 }
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.
double MiscGaussianDensity(double x, double mu, double sigma)
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 1612 of file ScalarSequence.C.

References queso_require_equal_to_msg, and queso_require_greater_equal_msg.

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

1618 {
1619  queso_require_equal_to_msg(centers.size(), bins.size(), "vectors 'centers' and 'bins' have different sizes");
1620 
1621  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1622 
1623  if (initialPos) {}; // just to remove compiler warning
1624 
1625  for (unsigned int j = 0; j < bins.size(); ++j) {
1626  centers[j] = 0.;
1627  bins[j] = 0;
1628  }
1629 
1630  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1631 
1632  double minCenter = minHorizontalValue - horizontalDelta/2.;
1633  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1634  for (unsigned int j = 0; j < centers.size(); ++j) {
1635  double factor = ((double) j)/(((double) centers.size()) - 1.);
1636  centers[j] = (1. - factor) * minCenter + factor * maxCenter;
1637  }
1638 
1639  unsigned int dataSize = this->subSequenceSize();
1640  for (unsigned int j = 0; j < dataSize; ++j) {
1641  double value = m_seq[j];
1642  if (value < minHorizontalValue) {
1643  bins[0]++;
1644  }
1645  else if (value >= maxHorizontalValue) {
1646  bins[bins.size()-1]++;
1647  }
1648  else {
1649  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1650  bins[index]++;
1651  }
1652  }
1653 
1654  return;
1655 }
std::vector< T > m_seq
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:90
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 2006 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().

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

1517 {
1518  queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid input");
1519 
1521  std::advance(pos1,initialPos);
1522 
1524  std::advance(pos2,initialPos+numPos);
1525 
1526  if ((initialPos+numPos) == this->subSequenceSize()) {
1527  queso_require_msg(!(pos2 != m_seq.end()), "invalid state");
1528  }
1529 
1531  pos = std::min_element(pos1, pos2);
1532  minValue = *pos;
1533  pos = std::max_element(pos1, pos2);
1534  maxValue = *pos;
1535 
1536  return;
1537 }
std::vector< T > m_seq
#define queso_require_less_equal_msg(expr1, expr2, msg)
Definition: asserts.h:89
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
std::vector< T >::const_iterator seqScalarPositionConstIteratorTypedef
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 }
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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...
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 2520 of file ScalarSequence.C.

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

2523 {
2524  queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
2525 
2526  T subMaxValue = subCorrespondingScalarValues.subMaxPlain();
2527  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
2528 
2529  unsigned int subNumPos = 0;
2530  for (unsigned int i = 0; i < iMax; ++i) {
2531  if (subCorrespondingScalarValues[i] == subMaxValue) {
2532  subNumPos++;
2533  }
2534  }
2535 
2536  subPositionsOfMaximum.resizeSequence(subNumPos);
2537  unsigned int j = 0;
2538  for (unsigned int i = 0; i < iMax; ++i) {
2539  if (subCorrespondingScalarValues[i] == subMaxValue) {
2540  subPositionsOfMaximum[j] = (*this)[i];
2541  j++;
2542  }
2543  }
2544 
2545  return subMaxValue;
2546 }
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:85
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.
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...
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 2210 of file ScalarSequence.C.

References queso_require_msg.

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

2214 {
2215  bool bRC = (initialPos < this->subSequenceSize());
2216  queso_require_msg(bRC, "invalid input data");
2217 
2218  unsigned int dataSize = this->subSequenceSize() - initialPos;
2219 
2220  T meanValue = this->subMeanExtra(initialPos,
2221  dataSize);
2222 
2223  T samValue = this->subSampleVarianceExtra(initialPos,
2224  dataSize,
2225  meanValue);
2226 
2227  T scaleValue;
2228  if (iqrValue <= 0.) {
2229  scaleValue = 1.06*std::sqrt(samValue)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
2230  }
2231  else {
2232  scaleValue = 1.06*std::min(std::sqrt(samValue),iqrValue/1.34)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
2233  }
2234 
2235  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2236  *m_env.subDisplayFile() << "In ScalarSequence<T>::subScaleForKde()"
2237  << ": iqrValue = " << iqrValue
2238  << ", meanValue = " << meanValue
2239  << ", samValue = " << samValue
2240  << ", dataSize = " << dataSize
2241  << ", scaleValue = " << scaleValue
2242  << std::endl;
2243  }
2244 
2245  return scaleValue;
2246 }
unsigned int displayVerbosity() const
Definition: Environment.C:400
const BaseEnvironment & m_env
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
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...
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 >
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 1885 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().

1888 {
1889  unsigned int numPos = this->subSequenceSize() - initialPos;
1890  sortedSequence.resizeSequence(numPos);
1891  this->extractScalarSeq(initialPos,
1892  1,
1893  numPos,
1894  sortedSequence);
1895  sortedSequence.subSort();
1896 
1897  return;
1898 }
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
void extractScalarSeq(unsigned int initialPos, unsigned int spacing, unsigned int numPos, ScalarSequence< T > &scalarSeq) const
Extracts a 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 3266 of file ScalarSequence.C.

3267 {
3268  std::sort(m_seq.begin(), m_seq.end());
3269  return;
3270 }
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 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.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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...
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 }
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
void 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.
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 }
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
void 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.
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 1791 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

1797 {
1798  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1799 
1800  for (unsigned int j = 0; j < bins.size(); ++j) {
1801  bins[j] = 0;
1802  }
1803 
1804  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1805  double minCenter = minHorizontalValue - horizontalDelta/2.;
1806  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1807  gridValues = new UniformOneDGrid<T>(m_env,
1808  "",
1809  bins.size(),
1810  minCenter,
1811  maxCenter);
1812 
1813  unsigned int dataSize = this->subSequenceSize();
1814  for (unsigned int j = 0; j < dataSize; ++j) {
1815  double value = m_seq[j];
1816  if (value < minHorizontalValue) {
1817  bins[0] += value;
1818  }
1819  else if (value >= maxHorizontalValue) {
1820  bins[bins.size()-1] += value;
1821  }
1822  else {
1823  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1824  bins[index] += value;
1825  }
1826  }
1827 
1828  return;
1829 }
std::vector< T > m_seq
const BaseEnvironment & m_env
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:90
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 1833 of file ScalarSequence.C.

References queso_require_greater_equal_msg.

1839 {
1840  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1841 
1842  for (unsigned int j = 0; j < bins.size(); ++j) {
1843  bins[j] = 0;
1844  }
1845 
1846  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1847  double minCenter = minHorizontalValue - horizontalDelta/2.;
1848  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1849  UniformOneDGrid<T> tmpGrid(m_env,
1850  "",
1851  bins.size(),
1852  minCenter,
1853  maxCenter);
1854  gridValues.clear();
1855  gridValues.resize(tmpGrid.size(),0.);
1856  for (unsigned int i = 0; i < tmpGrid.size(); ++i) {
1857  gridValues[i] = tmpGrid[i];
1858  }
1859 
1860  unsigned int dataSize = this->subSequenceSize();
1861  for (unsigned int j = 0; j < dataSize; ++j) {
1862  double value = m_seq[j];
1863  if (value < minHorizontalValue) {
1864  bins[0] += value;
1865  }
1866  else if (value >= maxHorizontalValue) {
1867  bins[bins.size()-1] += value;
1868  }
1869  else {
1870  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1871  bins[index] += value;
1872  }
1873  }
1874 
1875  //std::cout << "In ScalarSequence<T>::subWeightHistogram():" << std::endl;
1876  //for (unsigned int j = 0; j < bins.size(); ++j) {
1877  // std::cout << "bins[" << j << "] = " << bins[j] << std::endl;
1878  //}
1879 
1880  return;
1881 }
std::vector< T > m_seq
const BaseEnvironment & m_env
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:90
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 2580 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().

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

References UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, and UQ_FILE_EXTENSION_FOR_TXT_FORMAT.

2660 {
2661  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2662  this->writeSubMatlabHeader(ofs, this->subSequenceSize());
2663  }
2664  else if (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) {
2665  this->writeTxtHeader(ofs, this->subSequenceSize());
2666  }
2667 
2668  unsigned int chainSize = this->subSequenceSize();
2669  for (unsigned int j = 0; j < chainSize; ++j) {
2670  ofs << m_seq[j]
2671  << std::endl;
2672  }
2673 
2674  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2675  ofs << "];\n";
2676  }
2677 
2678  return;
2679 }
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:102
#define UQ_FILE_EXTENSION_FOR_TXT_FORMAT
Definition: Defines.h:103
std::vector< T > m_seq
void writeSubMatlabHeader(std::ofstream &ofs, double sequenceSize) const
Helper function to write header info for matlab files from one chain.
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 2351 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().

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

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

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

1547 {
1548  if (m_env.numSubEnvironments() == 1) {
1549  return this->subMinMaxExtra(initialPos,
1550  numPos,
1551  unifiedMinValue,
1552  unifiedMaxValue);
1553  }
1554 
1555  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1556 
1557  if (useOnlyInter0Comm) {
1558  if (m_env.inter0Rank() >= 0) {
1559  // Find local min and max
1560  T minValue;
1561  T maxValue;
1562  this->subMinMaxExtra(initialPos,
1563  numPos,
1564  minValue,
1565  maxValue);
1566 
1567  // Get overall min
1568  std::vector<double> sendBuf(1,0.);
1569  for (unsigned int i = 0; i < sendBuf.size(); ++i) {
1570  sendBuf[i] = minValue;
1571  }
1572  m_env.inter0Comm().template Allreduce<double>(&sendBuf[0], &unifiedMinValue, (int) sendBuf.size(), RawValue_MPI_MIN,
1573  "ScalarSequence<T>::unifiedMinMaxExtra()",
1574  "failed MPI.Allreduce() for min");
1575 
1576  // Get overall max
1577  for (unsigned int i = 0; i < sendBuf.size(); ++i) {
1578  sendBuf[i] = maxValue;
1579  }
1580  m_env.inter0Comm().template Allreduce<double>(&sendBuf[0], &unifiedMaxValue, (int) sendBuf.size(), RawValue_MPI_MAX,
1581  "ScalarSequence<T>::unifiedMinMaxExtra()",
1582  "failed MPI.Allreduce() for max");
1583 
1584  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1585  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMinMaxExtra()"
1586  << ": localMinValue = " << minValue
1587  << ", localMaxValue = " << maxValue
1588  << ", unifiedMinValue = " << unifiedMinValue
1589  << ", unifiedMaxValue = " << unifiedMaxValue
1590  << std::endl;
1591  }
1592  }
1593  else {
1594  // Node not in the 'inter0' communicator
1595  this->subMinMaxExtra(initialPos,
1596  numPos,
1597  unifiedMinValue,
1598  unifiedMaxValue);
1599  }
1600  }
1601  else {
1602  queso_error_msg("parallel vectors not supported yet");
1603  }
1604 
1605  //m_env.fullComm().Barrier();
1606 
1607  return;
1608 }
unsigned int displayVerbosity() const
Definition: Environment.C:400
#define queso_error_msg(msg)
Definition: asserts.h:47
const BaseEnvironment & m_env
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:292
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:265
#define RawValue_MPI_MIN
Definition: MpiComm.h:69
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
#define RawValue_MPI_MAX
Definition: MpiComm.h:70
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:271
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...
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 queso_error_msg(msg)
Definition: asserts.h:47
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
std::vector< T > m_seq
const BaseEnvironment & m_env
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:292
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:265
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:271
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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 ...
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 2550 of file ScalarSequence.C.

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

2553 {
2554  queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
2555 
2556  T maxValue = subCorrespondingScalarValues.subMaxPlain();
2557  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
2558 
2559  unsigned int numPos = 0;
2560  for (unsigned int i = 0; i < iMax; ++i) {
2561  if (subCorrespondingScalarValues[i] == maxValue) {
2562  numPos++;
2563  }
2564  }
2565 
2566  unifiedPositionsOfMaximum.resizeSequence(numPos);
2567  unsigned int j = 0;
2568  for (unsigned int i = 0; i < iMax; ++i) {
2569  if (subCorrespondingScalarValues[i] == maxValue) {
2570  unifiedPositionsOfMaximum[j] = (*this)[i];
2571  j++;
2572  }
2573  }
2574 
2575  return maxValue;
2576 }
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
T unifiedPositionsOfMaximum(const ScalarSequence< T > &subCorrespondingScalarValues, ScalarSequence< T > &unifiedPositionsOfMaximum)
Finds the positions where the maximum element occurs in the unified sequence.
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 2908 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().

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

References queso_error_msg, queso_require_msg, and RawValue_MPI_SUM.

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

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

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

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

2889 {
2890  ofs << m_name << "_sub" << m_env.subIdString() << " = zeros(" << sequenceSize
2891  << "," << 1
2892  << ");"
2893  << std::endl;
2894  ofs << m_name << "_sub" << m_env.subIdString() << " = [";
2895 }
const BaseEnvironment & m_env
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:305
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 2899 of file ScalarSequence.C.

2901 {
2902  ofs << sequenceSize << " " << 1 << std::endl;
2903 }
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 2875 of file ScalarSequence.C.

2877 {
2878  ofs << m_name << "_unified" << " = zeros(" << sequenceSize
2879  << "," << 1
2880  << ");"
2881  << std::endl;
2882  ofs << m_name << "_unified" << " = [";
2883 }

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 Fri Jun 17 2016 14:17:43 for queso-0.55.0 by  doxygen 1.8.5