25 #include <queso/MonteCarloSG.h> 
   26 #include <queso/GslVector.h> 
   27 #include <queso/GslMatrix.h> 
   32 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
   39   m_env                     (paramRv.env()),
 
   41   m_qoiFunction             (qoiFunction),
 
   42   m_paramSpace              (m_paramRv.imageSet().vectorSpace()),
 
   43   m_qoiSpace                (m_qoiFunction.imageSet().vectorSpace()),
 
   45   m_numPsNotSubWritten      (0),
 
   46   m_numQsNotSubWritten      (0),
 
   47   m_optionsObj              (alternativeOptionsValues),
 
   48   m_userDidNotProvideOptions(
false)
 
   50   if (m_env.subDisplayFile()) {
 
   51     *m_env.subDisplayFile() << 
"Entering MonteCarloSG<P_V,P_M,Q_V,Q_M>::constructor()" 
   52                             << 
": prefix = " << prefix
 
   53                             << 
", alternativeOptionsValues = " << alternativeOptionsValues
 
   54                             << 
", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
 
   59   if (m_optionsObj == NULL) {
 
   64     m_optionsObj = tempOptions;
 
   67     m_userDidNotProvideOptions = 
true;
 
   70   if (m_optionsObj->m_help != 
"") {
 
   71     if (m_env.subDisplayFile()) {
 
   72       *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
 
   76   queso_require_equal_to_msg(paramRv.imageSet().vectorSpace().dimLocal(), qoiFunction.domainSet().vectorSpace().dimLocal(), 
"'paramRv' and 'qoiFunction' are related to vector spaces of different dimensions");
 
   78   if (m_env.subDisplayFile()) {
 
   79     *m_env.subDisplayFile() << 
"Leaving MonteCarloSG<P_V,P_M,Q_V,Q_M>::constructor()" 
   84 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
   87   if (m_optionsObj && m_userDidNotProvideOptions) {
 
   91   if (m_qoiFunctionSynchronizer) 
delete m_qoiFunctionSynchronizer;
 
   94 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
  100   queso_require_equal_to_msg(m_qoiFunction.domainSet().vectorSpace().dimLocal(), workingPSeq.
vectorSizeLocal(), 
"'m_qoiFunction.domainSet' and 'workingPSeq' are related to vector spaces of different dimensions");
 
  102   queso_require_equal_to_msg(m_qoiFunction.imageSet().vectorSpace().dimLocal(), workingQSeq.
vectorSizeLocal(), 
"'m_qoiFunction.imageSet' and 'workingQSeq' are related to vector spaces of different dimensions");
 
  104   MiscCheckTheParallelEnvironment<P_V,Q_V>(m_paramRv.imageSet().vectorSpace().zeroVector(),
 
  105                                              m_qoiFunction.imageSet().vectorSpace().zeroVector());
 
  106   internGenerateSequence(m_paramRv,workingPSeq,workingQSeq);
 
  111 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
  118 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
  125   workingPSeq.
setName(m_optionsObj->m_prefix+
"ParamSeq");
 
  126   workingQSeq.
setName(m_optionsObj->m_prefix+
"QoiSeq");
 
  131   unsigned int subActualSizeBeforeGeneration = std::min(m_optionsObj->m_qseqSize,paramRv.
realizer().
subPeriod());
 
  132   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  133     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  134                             << 
": m_optionsObj->m_qseqSize = "                             << m_optionsObj->m_qseqSize
 
  136                             << 
", about to call actualGenerateSequence() with subActualSize = " << subActualSizeBeforeGeneration
 
  140     actualGenerateSequence(paramRv,
 
  143                            subActualSizeBeforeGeneration);
 
  146     actualReadSequence(paramRv,
 
  147                        m_optionsObj->m_qseqDataInputFileName,
 
  148                        m_optionsObj->m_qseqDataInputFileType,
 
  151                        subActualSizeBeforeGeneration);
 
  153   unsigned int subActualSizeAfterGeneration = workingPSeq.
subSequenceSize();
 
  156   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  157     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  158                             << 
": returned from call to actualGenerateSequence() with subActualSize = " << subActualSizeAfterGeneration
 
  165   if (m_env.subDisplayFile()) {
 
  166     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  167                             << 
", prefix = "                                                        << m_optionsObj->m_prefix
 
  168                             << 
": checking necessity of opening generic output file (qseq name is " << workingQSeq.
name()
 
  173   m_env.openOutputFile(m_optionsObj->m_dataOutputFileName,
 
  175                        m_optionsObj->m_dataOutputAllowedSet,
 
  184   if (m_env.subDisplayFile()) {
 
  185     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  186                             << 
", prefix = "                                            << m_optionsObj->m_prefix
 
  187                             << 
": checking necessity of opening output files for pseq " << workingPSeq.
name()
 
  193   if ((m_numPsNotSubWritten                        >  0                             ) &&
 
  195     workingPSeq.
subWriteContents(subActualSizeBeforeGeneration - m_numPsNotSubWritten,
 
  196                                  m_numPsNotSubWritten,
 
  197                                  m_optionsObj->m_pseqDataOutputFileName,
 
  198                                  m_optionsObj->m_pseqDataOutputFileType,
 
  199                                  m_optionsObj->m_pseqDataOutputAllowedSet);
 
  200     if (m_env.subDisplayFile()) {
 
  201       *m_env.subDisplayFile() << 
"In MonteCarloG<P_V,P_M>::internGenerateSequence()" 
  202                               << 
": just wrote remaining pseq positions (per period request)" 
  205     m_numPsNotSubWritten = 0;
 
  217     workingPSeq.
unifiedWriteContents(m_optionsObj->m_pseqDataOutputFileName,m_optionsObj->m_pseqDataOutputFileType);
 
  218     if (m_env.subDisplayFile()) {
 
  219       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  220                               << 
", prefix = "                         << m_optionsObj->m_prefix
 
  221                               << 
": closed unified data output file '" << m_optionsObj->m_pseqDataOutputFileName
 
  222                               << 
"' for pseq "                         << workingPSeq.
name()
 
  228 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS 
  229   if (m_optionsObj->m_pseqComputeStats) {
 
  230     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  231       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  232                               << 
": about to call 'workingPSeq.computeStatistics()'" 
  235     workingPSeq.computeStatistics(*m_optionsObj->m_pseqStatisticalOptionsObj,
 
  236                                   genericFilePtrSet.
ofsVar);
 
  237     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  238       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  239                               << 
": returned from call to 'workingPSeq.computeStatistics()'" 
  249   if (m_env.subDisplayFile()) {
 
  250     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  251                             << 
", prefix = "                                            << m_optionsObj->m_prefix
 
  252                             << 
": checking necessity of opening output files for qseq " << workingQSeq.
name()
 
  258   if ((m_numQsNotSubWritten                        >  0                             ) &&
 
  260     workingQSeq.
subWriteContents(subActualSizeBeforeGeneration - m_numQsNotSubWritten,
 
  261                                  m_numQsNotSubWritten,
 
  262                                  m_optionsObj->m_qseqDataOutputFileName,
 
  263                                  m_optionsObj->m_qseqDataOutputFileType,
 
  264                                  m_optionsObj->m_qseqDataOutputAllowedSet);
 
  265     if (m_env.subDisplayFile()) {
 
  266       *m_env.subDisplayFile() << 
"In MonteCarloG<P_V,P_M>::internGenerateSequence()" 
  267                               << 
": just wrote remaining qseq positions (per period request)" 
  270     m_numQsNotSubWritten = 0;
 
  282     workingQSeq.
unifiedWriteContents(m_optionsObj->m_qseqDataOutputFileName,m_optionsObj->m_qseqDataOutputFileType);
 
  283     if (m_env.subDisplayFile()) {
 
  284       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  285                               << 
", prefix = "                         << m_optionsObj->m_prefix
 
  286                               << 
": closed unified data output file '" << m_optionsObj->m_qseqDataOutputFileName
 
  287                               << 
"' for qseq "                         << workingQSeq.
name()
 
  293 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS 
  294   if (m_optionsObj->m_qseqComputeStats) {
 
  295     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  296       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  297                               << 
": about to call 'workingQSeq.computeStatistics()'" 
  300     workingQSeq.computeStatistics(*m_optionsObj->m_qseqStatisticalOptionsObj,
 
  301                                   genericFilePtrSet.
ofsVar);
 
  302     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  303       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  304                               << 
": returned from call to 'workingQSeq.computeStatistics()'" 
  312   if (genericFilePtrSet.
ofsVar) {
 
  314     delete genericFilePtrSet.
ofsVar;
 
  317     if (m_env.subDisplayFile()) {
 
  318       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  319                               << 
", prefix = "                         << m_optionsObj->m_prefix
 
  320                               << 
": closed generic data output file '" << m_optionsObj->m_dataOutputFileName
 
  321                               << 
"' for QoI sequence "                 << workingQSeq.
name()
 
  325   if (m_env.subDisplayFile()) {
 
  326     *m_env.subDisplayFile() << std::endl;
 
  332 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
  338         unsigned int                        requestedSeqSize)
 
  340   if (m_env.subDisplayFile()) {
 
  341     *m_env.subDisplayFile() << 
"Starting the generation of qoi sequence " << workingQSeq.
name()
 
  342                             << 
", with "                                  << requestedSeqSize
 
  348   struct timeval timevalSeq;
 
  349   struct timeval timevalQoIFunction;
 
  351   double seqRunTime         = 0;
 
  352   double qoiFunctionRunTime = 0;
 
  354   iRC = gettimeofday(&timevalSeq, NULL);
 
  358   m_numPsNotSubWritten = 0;
 
  360   m_numQsNotSubWritten = 0;
 
  362   P_V tmpP(m_paramSpace.zeroVector());
 
  363   Q_V tmpQ(m_qoiSpace.zeroVector());
 
  365   unsigned int actualSeqSize = 0;
 
  366   for (
unsigned int i = 0; i < requestedSeqSize; ++i) {
 
  369     if (m_optionsObj->m_qseqMeasureRunTimes) iRC = gettimeofday(&timevalQoIFunction, NULL);
 
  370     m_qoiFunctionSynchronizer->callFunction(&tmpP,NULL,&tmpQ,NULL,NULL,NULL); 
 
  371     if (m_optionsObj->m_qseqMeasureRunTimes) qoiFunctionRunTime += 
MiscGetEllapsedSeconds(&timevalQoIFunction);
 
  373     bool allQsAreFinite = 
true;
 
  374     for (
unsigned int j = 0; j < tmpQ.sizeLocal(); ++j) {
 
  375       if ((tmpQ[j] == INFINITY) || (tmpQ[j] == -INFINITY)) {
 
  376   std::cerr << 
"WARNING In MonteCarloSG<P_V,P_M,Q_V,Q_M>::actualGenerateSequence()" 
  377                   << 
", worldRank "      << m_env.worldRank()
 
  378                   << 
", fullRank "       << m_env.fullRank()
 
  379                   << 
", subEnvironment " << m_env.subId()
 
  380                   << 
", subRank "        << m_env.subRank()
 
  381                   << 
", inter0Rank "     << m_env.inter0Rank()
 
  383                   << 
", tmpQ[" << j << 
"] = " << tmpQ[j]
 
  384                   << 
", tmpP = "         << tmpP
 
  385                   << 
", tmpQ = "         << tmpQ
 
  387         allQsAreFinite = 
false;
 
  397     if (allQsAreFinite) {}; 
 
  401       m_numPsNotSubWritten++;
 
  402       if ((m_optionsObj->m_pseqDataOutputPeriod           >  0  ) &&
 
  403           (((i+1) % m_optionsObj->m_pseqDataOutputPeriod) == 0  ) &&
 
  404           (m_optionsObj->m_pseqDataOutputFileName         != 
".")) {
 
  406                                      m_optionsObj->m_pseqDataOutputPeriod,
 
  407                                      m_optionsObj->m_pseqDataOutputFileName,
 
  408                                      m_optionsObj->m_pseqDataOutputFileType,
 
  409                                      m_optionsObj->m_pseqDataOutputAllowedSet);
 
  410         if (m_env.subDisplayFile()) {
 
  411           *m_env.subDisplayFile() << 
"In MonteCarloG<P_V,P_M>::actualGenerateSequence()" 
  412                                   << 
": just wrote pseq positions (per period request)" 
  415         m_numPsNotSubWritten = 0;
 
  419       m_numQsNotSubWritten++;
 
  420       if ((m_optionsObj->m_qseqDataOutputPeriod           >  0  ) &&
 
  421           (((i+1) % m_optionsObj->m_qseqDataOutputPeriod) == 0  ) &&
 
  422           (m_optionsObj->m_qseqDataOutputFileName         != 
".")) {
 
  424                                      m_optionsObj->m_qseqDataOutputPeriod,
 
  425                                      m_optionsObj->m_qseqDataOutputFileName,
 
  426                                      m_optionsObj->m_qseqDataOutputFileType,
 
  427                                      m_optionsObj->m_qseqDataOutputAllowedSet);
 
  428         if (m_env.subDisplayFile()) {
 
  429           *m_env.subDisplayFile() << 
"In MonteCarloG<P_V,P_M>::actualGenerateSequence()" 
  430                                   << 
": just wrote qseq positions (per period request)" 
  433         m_numQsNotSubWritten = 0;
 
  440     if ((m_optionsObj->m_qseqDisplayPeriod            > 0) &&
 
  441         (((i+1) % m_optionsObj->m_qseqDisplayPeriod) == 0)) {
 
  442       if (m_env.subDisplayFile()) {
 
  443         *m_env.subDisplayFile() << 
"Finished generating " << i+1
 
  457   if (m_env.subDisplayFile()) {
 
  458     *m_env.subDisplayFile() << 
"Finished the generation of qoi sequence " << workingQSeq.
name()
 
  461                             << 
"\nSome information about this sequence:" 
  462                             << 
"\n  Sequence run time = " << seqRunTime
 
  464                             << 
"\n\n Breaking of the seq run time:\n" 
  465                             << 
"\n  QoI function run time   = " << qoiFunctionRunTime
 
  466                             << 
" seconds ("                     << 100.*qoiFunctionRunTime/seqRunTime
 
  474 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
  478   const std::string&                        dataInputFileName,
 
  479   const std::string&                        dataInputFileType,
 
  482         unsigned int                        requestedSeqSize)
 
  485   P_V tmpP(m_paramSpace.zeroVector());
 
  486   for (
unsigned int i = 0; i < requestedSeqSize; ++i) {
 
MonteCarloSG(const char *prefix, const McOptionsValues *alternativeOptionsValues, const BaseVectorRV< P_V, P_M > ¶mRv, const BaseVectorFunction< P_V, P_M, Q_V, Q_M > &qoiFunction)
Constructor. 
 
void actualGenerateSequence(const BaseVectorRV< P_V, P_M > ¶mRv, BaseVectorSequence< P_V, P_M > &workingPSeq, BaseVectorSequence< Q_V, Q_M > &workingQSeq, unsigned int seqSize)
This method actually generates the QoI sequence. 
 
void setName(const std::string &newName)
Changes the name of the sequence of vectors. 
 
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization. 
 
double MiscGetEllapsedSeconds(struct timeval *timeval0)
 
A templated class for synchronizing the calls of vector-valued functions. 
 
void generateSequence(BaseVectorSequence< P_V, P_M > &workingPSeq, BaseVectorSequence< Q_V, Q_M > &workingQSeq)
Generates the QoI (output) sequence, it calls internGenerateSequence(). 
 
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization. 
 
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization. 
 
unsigned int vectorSizeLocal() const 
Local dimension (size) of the vector space. 
 
A templated class that implements a Monte Carlo generator of samples. 
 
#define queso_require_equal_to_msg(expr1, expr2, msg)
 
virtual void subWriteContents(unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const =0
Writes info of the sub-sequence to a file. See template specialization. 
 
Struct for handling data input and output from files. 
 
virtual void realization(V &nextValues) const =0
Performs a realization (sample) from a probability density function. See template specialization...
 
void print(std::ostream &os) const 
Prints the sequence. 
 
~MonteCarloSG()
Destructor. 
 
virtual void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const =0
Writes info of the unified sequence to a file. See template specialization. 
 
void internGenerateSequence(const BaseVectorRV< P_V, P_M > ¶mRv, BaseVectorSequence< P_V, P_M > &workingPSeq, BaseVectorSequence< Q_V, Q_M > &workingQSeq)
Generates the QoI (output) sequence; it calls actualGenerateSequence(). 
 
std::ofstream * ofsVar
Provides a stream interface to write data to files. 
 
virtual void getPositionValues(unsigned int posId, V &vec) const =0
Gets the values of the sequence at position posId and stores them at vec. See template specialization...
 
virtual void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)=0
Reads info of the unified sequence from a file. See template specialization. 
 
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
 
void actualReadSequence(const BaseVectorRV< P_V, P_M > ¶mRv, const std::string &dataInputFileName, const std::string &dataInputFileType, BaseVectorSequence< P_V, P_M > &workingPSeq, BaseVectorSequence< Q_V, Q_M > &workingQSeq, unsigned int seqSize)
Reads the sequence. 
 
unsigned int subPeriod() const 
Sub-period of the realization. Access to protected attribute m_subPeriod. 
 
const BaseVectorRealizer< V, M > & realizer() const 
Finds a realization (sample) of the PDF of this vector RV; access to private attribute m_realizer...
 
#define UQ_MOC_SG_FILENAME_FOR_NO_FILE
 
const std::string & name() const 
Access to protected attribute m_name: name of the sequence of vectors. 
 
This class provides options for the Monte Carlo sequence generator if no input file is available...