MetropolisHastingsSG.h

Go to the documentation of this file.

//-----------------------------------------------------------------------bl-
//--------------------------------------------------------------------------
//
// QUESO - a library to support the Quantification of Uncertainty
// for Estimation, Simulation and Optimization
//
// Copyright (C) 2008,2009,2010,2011,2012,2013 The PECOS Development Team
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the Version 2.1 GNU Lesser General
// Public License as published by the Free Software Foundation.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc. 51 Franklin Street, Fifth Floor,
// Boston, MA  02110-1301  USA
//
//-----------------------------------------------------------------------el-

#ifndef UQ_MH_SG_H
#define UQ_MH_SG_H

#include <queso/MetropolisHastingsSGOptions.h>
#include <queso/TKGroup.h>
#include <queso/VectorRV.h>
#include <queso/VectorSpace.h>
#include <queso/MarkovChainPositionData.h>
#include <queso/ScalarFunctionSynchronizer.h>
#include <queso/SequenceOfVectors.h>
#include <queso/ArrayOfSequences.h>
#include <sys/time.h>
#include <fstream>
#include <boost/math/special_functions.hpp> // for Boost isnan. Note parentheses are important in function call.

namespace QUESO {

//--------------------------------------------------
// MHRawChainInfoStruct --------------------------
//--------------------------------------------------
struct MHRawChainInfoStruct
{

  MHRawChainInfoStruct();

  MHRawChainInfoStruct(const MHRawChainInfoStruct& rhs);

  ~MHRawChainInfoStruct();


  MHRawChainInfoStruct& operator= (const MHRawChainInfoStruct& rhs);

  MHRawChainInfoStruct& operator+=(const MHRawChainInfoStruct& rhs);


  void copy  (const MHRawChainInfoStruct& src);

  void reset ();

  void mpiSum(const MpiComm& comm, MHRawChainInfoStruct& sumInfo) const;

  double       runTime;
  double       candidateRunTime;
  double       targetRunTime;
  double       mhAlphaRunTime;
  double       drAlphaRunTime;
  double       drRunTime;
  double       amRunTime;

  unsigned int numTargetCalls;
  unsigned int numDRs;
  unsigned int numOutOfTargetSupport;
  unsigned int numOutOfTargetSupportInDR;
  unsigned int numRejections;

};

//--------------------------------------------------
// MetropolisHastingsSG----------------------
//--------------------------------------------------

template <class P_V,class P_M>
class MetropolisHastingsSG
{
public:


  MetropolisHastingsSG(const char*                         prefix,
                              const MhOptionsValues*       alternativeOptionsValues, // dakota
                              const BaseVectorRV<P_V,P_M>& sourceRv,
                              const P_V&                          initialPosition,
                              const P_M*                          inputProposalCovMatrix);

  MetropolisHastingsSG(const MLSamplingLevelOptions& mlOptions,
                              const BaseVectorRV<P_V,P_M>&  sourceRv,
                              const P_V&                           initialPosition,
                              const P_M*                           inputProposalCovMatrix);

  ~MetropolisHastingsSG();



  void         generateSequence   (BaseVectorSequence<P_V,P_M>& workingChain,
                                   ScalarSequence<double>*      workingLogLikelihoodValues,
                                   ScalarSequence<double>*      workingLogTargetValues);

  void         getRawChainInfo    (MHRawChainInfoStruct& info) const;




  void   print                    (std::ostream& os) const;
  friend std::ostream& operator<<(std::ostream& os,
      const MetropolisHastingsSG<P_V,P_M>& obj)
  {
    obj.print(os);

    return os;
  }

private:

  void   commonConstructor        ();


  void   generateFullChain        (const P_V&                          valuesOf1stPosition,
                                   unsigned int                        chainSize,
                                   BaseVectorSequence<P_V,P_M>& workingChain,
                                   ScalarSequence<double>*      workingLogLikelihoodValues,
                                   ScalarSequence<double>*      workingLogTargetValues);

  void   readFullChain            (const std::string&                  inputFileName,
                                   const std::string&                  inputFileType,
                                   unsigned int                        chainSize,
                                   BaseVectorSequence<P_V,P_M>& workingChain);


  void   updateAdaptedCovMatrix   (const BaseVectorSequence<P_V,P_M>&  subChain,
                                   unsigned int                               idOfFirstPositionInSubChain,
                                   double&                                    lastChainSize,
                                   P_V&                                       lastMean,
                                   P_M&                                       lastAdaptedCovMatrix);


  double alpha                    (const MarkovChainPositionData<P_V>& x,
                                   const MarkovChainPositionData<P_V>& y,
                                   unsigned int                               xStageId,
                                   unsigned int                               yStageId,
                                   double*                                    alphaQuotientPtr = NULL);


  double alpha                    (const std::vector<MarkovChainPositionData<P_V>*>& inputPositions,
                                   const std::vector<unsigned int                        >& inputTKStageIds);


  bool   acceptAlpha              (double                                     alpha);


  int    writeInfo                (const BaseVectorSequence<P_V,P_M>&  workingChain,
                                   std::ofstream&                             ofsvar) const;

  const BaseEnvironment&                     m_env;
  const VectorSpace <P_V,P_M>&               m_vectorSpace;
  const BaseJointPdf<P_V,P_M>&               m_targetPdf;
        P_V                                         m_initialPosition;
        P_M                                         m_initialProposalCovMatrix;
        bool                                        m_nullInputProposalCovMatrix;
        unsigned int                                m_numDisabledParameters; // gpmsa2
        std::vector<bool>                           m_parameterEnabledStatus; // gpmsa2
  const ScalarFunctionSynchronizer<P_V,P_M>* m_targetPdfSynchronizer;

        BaseTKGroup<P_V,P_M>*                m_tk;
        unsigned int                                m_positionIdForDebugging;
        unsigned int                                m_stageIdForDebugging;
        std::vector<unsigned int>                   m_idsOfUniquePositions;
        std::vector<double>                         m_logTargets;
        std::vector<double>                         m_alphaQuotients;
        double                                      m_lastChainSize;
        P_V*                                        m_lastMean;
        P_M*                                        m_lastAdaptedCovMatrix;
        unsigned int                                m_numPositionsNotSubWritten;

        MHRawChainInfoStruct                      m_rawChainInfo;

        MhOptionsValues                      m_alternativeOptionsValues;
        MetropolisHastingsSGOptions*         m_optionsObj;
};

}  // End namespace QUESO

#endif // UQ_MH_SG_H