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-2015 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 {

class GslVector;
class GslMatrix;

//--------------------------------------------------
// 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 = GslVector, class P_M = GslMatrix>
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 char*                         prefix,
                       const MhOptionsValues*       alternativeOptionsValues, // dakota
                       const BaseVectorRV<P_V,P_M>& sourceRv,
                       const P_V&                          initialPosition,
                       double                              initialLogPrior,
                       double                              initialLogLikelihood,
                       const P_M*                          inputProposalCovMatrix);

  MetropolisHastingsSG(const MLSamplingLevelOptions& mlOptions,
                       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,
                       double                               initialLogPrior,
                       double                               initialLogLikelihood,
                       const P_M*                           inputProposalCovMatrix);

  ~MetropolisHastingsSG();



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

  void         getRawChainInfo    (MHRawChainInfoStruct& info) const;


  const BaseTKGroup<P_V, P_M> & transitionKernel() 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 adapt(unsigned int positionId,
      BaseVectorSequence<P_V, P_M> & workingChain);

  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;

  const MhOptionsValues * m_optionsObj;
  MetropolisHastingsSGOptions * m_oldOptions;

        bool m_computeInitialPriorAndLikelihoodValues;
        double m_initialLogPriorValue;
        double m_initialLogLikelihoodValue;

  void transformInitialCovMatrixToGaussianSpace(const BoxSubset<P_V, P_M> &
      boxSubset);

  bool m_userDidNotProvideOptions;
};

}  // End namespace QUESO

#endif // UQ_MH_SG_H