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 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS 
   48   m_alternativeOptionsValues(NULL,NULL),
 
   50   m_alternativeOptionsValues(),
 
   54   if (m_env.subDisplayFile()) {
 
   55     *m_env.subDisplayFile() << 
"Entering MonteCarloSG<P_V,P_M,Q_V,Q_M>::constructor()" 
   56                             << 
": prefix = " << prefix
 
   57                             << 
", alternativeOptionsValues = " << alternativeOptionsValues
 
   58                             << 
", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
 
   62   if (alternativeOptionsValues) m_alternativeOptionsValues = *alternativeOptionsValues;
 
   63   if (m_env.optionsInputFileName() == 
"") {
 
   68     m_optionsObj->scanOptionsValues();
 
   71   UQ_FATAL_TEST_MACRO(paramRv.imageSet().vectorSpace().dimLocal() != qoiFunction.domainSet().vectorSpace().dimLocal(),
 
   73                       "MonteCarloSG<P_V,P_M,Q_V,Q_M>::constructor()",
 
   74                       "'paramRv' and 'qoiFunction' are related to vector spaces of different dimensions");
 
   76   if (m_env.subDisplayFile()) {
 
   77     *m_env.subDisplayFile() << 
"Leaving MonteCarloSG<P_V,P_M,Q_V,Q_M>::constructor()" 
   82 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
   85   if (m_optionsObj             ) 
delete m_optionsObj;
 
   86   if (m_qoiFunctionSynchronizer) 
delete m_qoiFunctionSynchronizer;
 
   89 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
   97                       "MonteCarloSG<P_V,P_M,Q_V,Q_M>::generateSequence()",
 
   98                       "'m_qoiFunction.domainSet' and 'workingPSeq' are related to vector spaces of different dimensions");
 
  102                       "MonteCarloSG<P_V,P_M,Q_V,Q_M>::generateSequence()",
 
  103                       "'m_qoiFunction.imageSet' and 'workingQSeq' are related to vector spaces of different dimensions");
 
  105   MiscCheckTheParallelEnvironment<P_V,Q_V>(m_paramRv.imageSet().vectorSpace().zeroVector(),
 
  106                                              m_qoiFunction.imageSet().vectorSpace().zeroVector());
 
  107   internGenerateSequence(m_paramRv,workingPSeq,workingQSeq);
 
  112 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
  119 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
  126   workingPSeq.
setName(m_optionsObj->m_prefix+
"ParamSeq");
 
  127   workingQSeq.
setName(m_optionsObj->m_prefix+
"QoiSeq");
 
  132   unsigned int subActualSizeBeforeGeneration = std::min(m_optionsObj->m_ov.m_qseqSize,paramRv.
realizer().
subPeriod());
 
  133   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  134     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  135                             << 
": m_optionsObj->m_ov.m_qseqSize = "                             << m_optionsObj->m_ov.m_qseqSize
 
  137                             << 
", about to call actualGenerateSequence() with subActualSize = " << subActualSizeBeforeGeneration
 
  141     actualGenerateSequence(paramRv,
 
  144                            subActualSizeBeforeGeneration);
 
  147     actualReadSequence(paramRv,
 
  148                        m_optionsObj->m_ov.m_qseqDataInputFileName,
 
  149                        m_optionsObj->m_ov.m_qseqDataInputFileType,
 
  152                        subActualSizeBeforeGeneration);
 
  154   unsigned int subActualSizeAfterGeneration = workingPSeq.
subSequenceSize();
 
  157                       "MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()",
 
  158                       "P and Q sequences should have the same size!");
 
  160   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  161     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  162                             << 
": returned from call to actualGenerateSequence() with subActualSize = " << subActualSizeAfterGeneration
 
  169   if (m_env.subDisplayFile()) {
 
  170     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  171                             << 
", prefix = "                                                        << m_optionsObj->m_prefix
 
  172                             << 
": checking necessity of opening generic output file (qseq name is " << workingQSeq.
name()
 
  177   m_env.openOutputFile(m_optionsObj->m_ov.m_dataOutputFileName,
 
  179                        m_optionsObj->m_ov.m_dataOutputAllowedSet,
 
  188   if (m_env.subDisplayFile()) {
 
  189     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  190                             << 
", prefix = "                                            << m_optionsObj->m_prefix
 
  191                             << 
": checking necessity of opening output files for pseq " << workingPSeq.
name()
 
  197   if ((m_numPsNotSubWritten                        >  0                             ) &&
 
  199     workingPSeq.
subWriteContents(subActualSizeBeforeGeneration - m_numPsNotSubWritten,
 
  200                                  m_numPsNotSubWritten,
 
  201                                  m_optionsObj->m_ov.m_pseqDataOutputFileName,
 
  202                                  m_optionsObj->m_ov.m_pseqDataOutputFileType,
 
  203                                  m_optionsObj->m_ov.m_pseqDataOutputAllowedSet);
 
  204     if (m_env.subDisplayFile()) {
 
  205       *m_env.subDisplayFile() << 
"In MonteCarloG<P_V,P_M>::internGenerateSequence()" 
  206                               << 
": just wrote remaining pseq positions (per period request)" 
  209     m_numPsNotSubWritten = 0;
 
  221     workingPSeq.
unifiedWriteContents(m_optionsObj->m_ov.m_pseqDataOutputFileName,m_optionsObj->m_ov.m_pseqDataOutputFileType);
 
  222     if (m_env.subDisplayFile()) {
 
  223       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  224                               << 
", prefix = "                         << m_optionsObj->m_prefix
 
  225                               << 
": closed unified data output file '" << m_optionsObj->m_ov.m_pseqDataOutputFileName
 
  226                               << 
"' for pseq "                         << workingPSeq.
name()
 
  232 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS 
  233   if (m_optionsObj->m_ov.m_pseqComputeStats) {
 
  234     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  235       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  236                               << 
": about to call 'workingPSeq.computeStatistics()'" 
  239     workingPSeq.computeStatistics(*m_optionsObj->m_pseqStatisticalOptionsObj,
 
  240                                   genericFilePtrSet.
ofsVar);
 
  241     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  242       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  243                               << 
": returned from call to 'workingPSeq.computeStatistics()'" 
  253   if (m_env.subDisplayFile()) {
 
  254     *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  255                             << 
", prefix = "                                            << m_optionsObj->m_prefix
 
  256                             << 
": checking necessity of opening output files for qseq " << workingQSeq.
name()
 
  262   if ((m_numQsNotSubWritten                        >  0                             ) &&
 
  264     workingQSeq.
subWriteContents(subActualSizeBeforeGeneration - m_numQsNotSubWritten,
 
  265                                  m_numQsNotSubWritten,
 
  266                                  m_optionsObj->m_ov.m_qseqDataOutputFileName,
 
  267                                  m_optionsObj->m_ov.m_qseqDataOutputFileType,
 
  268                                  m_optionsObj->m_ov.m_qseqDataOutputAllowedSet);
 
  269     if (m_env.subDisplayFile()) {
 
  270       *m_env.subDisplayFile() << 
"In MonteCarloG<P_V,P_M>::internGenerateSequence()" 
  271                               << 
": just wrote remaining qseq positions (per period request)" 
  274     m_numQsNotSubWritten = 0;
 
  286     workingQSeq.
unifiedWriteContents(m_optionsObj->m_ov.m_qseqDataOutputFileName,m_optionsObj->m_ov.m_qseqDataOutputFileType);
 
  287     if (m_env.subDisplayFile()) {
 
  288       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  289                               << 
", prefix = "                         << m_optionsObj->m_prefix
 
  290                               << 
": closed unified data output file '" << m_optionsObj->m_ov.m_qseqDataOutputFileName
 
  291                               << 
"' for qseq "                         << workingQSeq.
name()
 
  297 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS 
  298   if (m_optionsObj->m_ov.m_qseqComputeStats) {
 
  299     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  300       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  301                               << 
": about to call 'workingQSeq.computeStatistics()'" 
  304     workingQSeq.computeStatistics(*m_optionsObj->m_qseqStatisticalOptionsObj,
 
  305                                   genericFilePtrSet.
ofsVar);
 
  306     if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
 
  307       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  308                               << 
": returned from call to 'workingQSeq.computeStatistics()'" 
  316   if (genericFilePtrSet.
ofsVar) {
 
  318     delete genericFilePtrSet.
ofsVar;
 
  321     if (m_env.subDisplayFile()) {
 
  322       *m_env.subDisplayFile() << 
"In MonteCarloSG<P_V,P_M,Q_V,Q_M>::internGenerateSequence()" 
  323                               << 
", prefix = "                         << m_optionsObj->m_prefix
 
  324                               << 
": closed generic data output file '" << m_optionsObj->m_ov.m_dataOutputFileName
 
  325                               << 
"' for QoI sequence "                 << workingQSeq.
name()
 
  329   if (m_env.subDisplayFile()) {
 
  330     *m_env.subDisplayFile() << std::endl;
 
  336 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
  342         unsigned int                        requestedSeqSize)
 
  344   if (m_env.subDisplayFile()) {
 
  345     *m_env.subDisplayFile() << 
"Starting the generation of qoi sequence " << workingQSeq.
name()
 
  346                             << 
", with "                                  << requestedSeqSize
 
  352   struct timeval timevalSeq;
 
  353   struct timeval timevalQoIFunction;
 
  355   double seqRunTime         = 0;
 
  356   double qoiFunctionRunTime = 0;
 
  358   iRC = gettimeofday(&timevalSeq, NULL);
 
  362   m_numPsNotSubWritten = 0;
 
  364   m_numQsNotSubWritten = 0;
 
  366   P_V tmpP(m_paramSpace.zeroVector());
 
  367   Q_V tmpQ(m_qoiSpace.zeroVector());
 
  369   unsigned int actualSeqSize = 0;
 
  370   for (
unsigned int i = 0; i < requestedSeqSize; ++i) {
 
  373     if (m_optionsObj->m_ov.m_qseqMeasureRunTimes) iRC = gettimeofday(&timevalQoIFunction, NULL);
 
  374     m_qoiFunctionSynchronizer->callFunction(&tmpP,NULL,&tmpQ,NULL,NULL,NULL); 
 
  375     if (m_optionsObj->m_ov.m_qseqMeasureRunTimes) qoiFunctionRunTime += 
MiscGetEllapsedSeconds(&timevalQoIFunction);
 
  377     bool allQsAreFinite = 
true;
 
  378     for (
unsigned int j = 0; j < tmpQ.sizeLocal(); ++j) {
 
  379       if ((tmpQ[j] == INFINITY) || (tmpQ[j] == -INFINITY)) {
 
  380   std::cerr << 
"WARNING In MonteCarloSG<P_V,P_M,Q_V,Q_M>::actualGenerateSequence()" 
  381                   << 
", worldRank "      << m_env.worldRank()
 
  382                   << 
", fullRank "       << m_env.fullRank()
 
  383                   << 
", subEnvironment " << m_env.subId()
 
  384                   << 
", subRank "        << m_env.subRank()
 
  385                   << 
", inter0Rank "     << m_env.inter0Rank()
 
  387                   << 
", tmpQ[" << j << 
"] = " << tmpQ[j]
 
  388                   << 
", tmpP = "         << tmpP
 
  389                   << 
", tmpQ = "         << tmpQ
 
  391         allQsAreFinite = 
false;
 
  401     if (allQsAreFinite) {}; 
 
  405       m_numPsNotSubWritten++;
 
  406       if ((m_optionsObj->m_ov.m_pseqDataOutputPeriod           >  0  ) &&
 
  407           (((i+1) % m_optionsObj->m_ov.m_pseqDataOutputPeriod) == 0  ) &&
 
  408           (m_optionsObj->m_ov.m_pseqDataOutputFileName         != 
".")) {
 
  409         workingPSeq.
subWriteContents(i + 1 - m_optionsObj->m_ov.m_pseqDataOutputPeriod,
 
  410                                      m_optionsObj->m_ov.m_pseqDataOutputPeriod,
 
  411                                      m_optionsObj->m_ov.m_pseqDataOutputFileName,
 
  412                                      m_optionsObj->m_ov.m_pseqDataOutputFileType,
 
  413                                      m_optionsObj->m_ov.m_pseqDataOutputAllowedSet);
 
  414         if (m_env.subDisplayFile()) {
 
  415           *m_env.subDisplayFile() << 
"In MonteCarloG<P_V,P_M>::actualGenerateSequence()" 
  416                                   << 
": just wrote pseq positions (per period request)" 
  419         m_numPsNotSubWritten = 0;
 
  423       m_numQsNotSubWritten++;
 
  424       if ((m_optionsObj->m_ov.m_qseqDataOutputPeriod           >  0  ) &&
 
  425           (((i+1) % m_optionsObj->m_ov.m_qseqDataOutputPeriod) == 0  ) &&
 
  426           (m_optionsObj->m_ov.m_qseqDataOutputFileName         != 
".")) {
 
  427         workingQSeq.
subWriteContents(i + 1 - m_optionsObj->m_ov.m_qseqDataOutputPeriod,
 
  428                                      m_optionsObj->m_ov.m_qseqDataOutputPeriod,
 
  429                                      m_optionsObj->m_ov.m_qseqDataOutputFileName,
 
  430                                      m_optionsObj->m_ov.m_qseqDataOutputFileType,
 
  431                                      m_optionsObj->m_ov.m_qseqDataOutputAllowedSet);
 
  432         if (m_env.subDisplayFile()) {
 
  433           *m_env.subDisplayFile() << 
"In MonteCarloG<P_V,P_M>::actualGenerateSequence()" 
  434                                   << 
": just wrote qseq positions (per period request)" 
  437         m_numQsNotSubWritten = 0;
 
  444     if ((m_optionsObj->m_ov.m_qseqDisplayPeriod            > 0) &&
 
  445         (((i+1) % m_optionsObj->m_ov.m_qseqDisplayPeriod) == 0)) {
 
  446       if (m_env.subDisplayFile()) {
 
  447         *m_env.subDisplayFile() << 
"Finished generating " << i+1
 
  461   if (m_env.subDisplayFile()) {
 
  462     *m_env.subDisplayFile() << 
"Finished the generation of qoi sequence " << workingQSeq.
name()
 
  465                             << 
"\nSome information about this sequence:" 
  466                             << 
"\n  Sequence run time = " << seqRunTime
 
  468                             << 
"\n\n Breaking of the seq run time:\n" 
  469                             << 
"\n  QoI function run time   = " << qoiFunctionRunTime
 
  470                             << 
" seconds ("                     << 100.*qoiFunctionRunTime/seqRunTime
 
  478 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
 
  482   const std::string&                        dataInputFileName,
 
  483   const std::string&                        dataInputFileType,
 
  486         unsigned int                        requestedSeqSize)
 
  489   P_V tmpP(m_paramSpace.zeroVector());
 
  490   for (
unsigned int i = 0; i < requestedSeqSize; ++i) {
 
void setName(const std::string &newName)
Changes the name of the sequence of vectors. 
double MiscGetEllapsedSeconds(struct timeval *timeval0)
#define UQ_MOC_SG_FILENAME_FOR_NO_FILE
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization. 
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. 
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. 
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. 
const BaseVectorRealizer< V, M > & realizer() const 
Finds a realization (sample) of the PDF of this vector RV; access to private attribute m_realizer...
A templated class that implements a Monte Carlo generator of samples. 
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 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...
This class provides options for the Monte Carlo sequence generator if no input file is available...
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 vectorSizeLocal() const 
Local dimension (size) of the vector space. 
std::ofstream * ofsVar
Provides a stream interface to write data to files. 
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
~MonteCarloSG()
Destructor. 
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(). 
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. 
A templated class for synchronizing the calls of vector-valued functions. 
const std::string & name() const 
Access to protected attribute m_name: name of the sequence of vectors. 
virtual void realization(V &nextValues) const =0
Performs a realization (sample) from a probability density function. See template specialization...
This class reads the options for the Monte Carlo sequence generator from an input file...
void print(std::ostream &os) const 
Prints the sequence. 
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization. 
void generateSequence(BaseVectorSequence< P_V, P_M > &workingPSeq, BaseVectorSequence< Q_V, Q_M > &workingQSeq)
Generates the QoI (output) sequence, it calls internGenerateSequence(). 
unsigned int subPeriod() const 
Sub-period of the realization. Access to protected attribute m_subPeriod. 
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.