mcmc.h

Go to the documentation of this file.

/* =====================================================================================
 
                      The UQ Toolkit (UQTk) version 3.1.5
                          Copyright (2024) NTESS
                        https://www.sandia.gov/UQToolkit/
                        https://github.com/sandialabs/UQTk
 
     Copyright 2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
     Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government
     retains certain rights in this software.
 
     This file is part of The UQ Toolkit (UQTk)
 
     UQTk is open source software: you can redistribute it and/or modify
     it under the terms of BSD 3-Clause License
 
     UQTk 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
     BSD 3 Clause License for more details.
 
     You should have received a copy of the BSD 3 Clause License
     along with UQTk. If not, see https://choosealicense.com/licenses/bsd-3-clause/.
 
     Questions? Contact the UQTk Developers at https://github.com/sandialabs/UQTk/discussions
     Sandia National Laboratories, Livermore, CA, USA
===================================================================================== */
 
#ifndef UQTKMCMC_H_SEEN
#define UQTKMCMC_H_SEEN
 
#include "dsfmt_add.h"
#include "arrayio.h"
#include "arraytools.h"
 
#include <iostream>
#include <string.h>
#include <stdio.h>
#include <sstream>
 
using namespace std; // needed for python string conversion
 
class LogPosteriorBase{
public:
  virtual double eval(Array1D<double>&){
      return 3.14;
  };
  virtual ~LogPosteriorBase(){};
};
 
class base{
public:
    virtual double fun(double* x, int n){return 0.0;}
};
 
class main{
public:
    base* b_;
 
    main(base& b){
        b_ = &b;
    }
};
//*****************************************
 
class MCMC{
public:
    // Struct for the chain state
    struct chainstate{
      int step;
      Array1D<double> state;
      double alfa;
      double post;
    };
 
    // Constructors:
 
    MCMC(double (*logposterior)(Array1D<double>&, void *), void *postinfo);
    MCMC(LogPosteriorBase& L);
    MCMC();
 
    // Set or initialization functions:
 
    void setWriteFlag(int I);
    void setFcnAccept(void (*fcnAccept)(void *));
    void setFcnReject(void (*fcnReject)(void *));
    void setChainDim(int chdim);
    void initChainPropCov(Array2D<double>& propcov);
    void initChainPropCovDiag(Array1D<double>& sig);
    void setOutputInfo(string outtype, string file,int freq_file, int freq_screen);
    void namesPrepended();
    void setSeed(int seed);
    void setLower(double lower, int i);
    void setUpper(double upper, int i);
    void setDefaultDomain();
    void setPostInfo(void *postinfo);
 
    // Get functions:
 
    void getChainPropCov(Array2D<double>& propcov);
    string getFilename();
    int getWriteFlag();
    void getSamples(int burnin, int every,Array2D<double>& samples);
    void getSamples(Array2D<double>& samples);
    void getFcnAccept(void (*fcnAccept)(void *));
    void getFcnReject(void (*fcnReject)(void *));
    string getOutputType();
    int getFileFreq();
    int getScreenFreq();
    bool getNamesPrepended();
    int getSeed();
    double getLower(int i);
    double getUpper(int i);
    bool getDimInit();
    void getPostInfo(void *post);
    bool getPropCovInit();
    bool getOutputInit();
    int getLastWrite();
    bool getFcnAcceptInit();
    bool getFcnRejectInit();
    virtual int getNSubSteps(){return 1;};
    int getLowerFlag(int i);
    int getUpperFlag(int i);
    void getAcceptRatio(double * accrat);
    double getAcceptRatio();
    int GetChainDim() const;
 
    // Chain Functions:
 
    void resetChainState();
    void resetChainFilename(string filename);
    void parseBinChain(string filename, Array1D<chainstate>& readchain);
    void writeFullChainTxt(string filename, Array1D<chainstate> fullchain);
    void getFullChain(Array1D<chainstate>& readchain);
    void appendMAP();
    double getMode(Array1D<double>& MAPparams);
    int getFullChainSize();
 
    // Functions to make sure the code respects the interface for chainstates:
 
    void setCurrentStateStep(int i);
    void getCurrentStateState(Array1D<double>& state);
    double getCurrentStatePost();
    void setCurrentStateState(Array1D<double>& newState);
    void setCurrentStatePost(double newPost);
    void setCurrentStateAlfa(double newAlfa);
    double getModeStatePost();
    void getModeStateState(Array1D<double>& state);
 
    // Run functions:
 
    virtual void runOptim(Array1D<double>& start);
    virtual void runChain(int ncalls, Array1D<double>& chstart) = 0;
    virtual void runChain(int ncalls) = 0;
    void runAcceptFcn();
    void runRejectFcn();
 
 
 
    bool newModeFound();
    double evalLogPosterior(Array1D<double>& m);
    bool inDomain(Array1D<double>& m);
    void writeChainTxt(string filename);
    void writeChainBin(string filename);
 
    void setNewMode(bool value);
 
 
 
protected:
    void setAcceptRatio(double d);
    void addCurrentState();
    void updateMode();
    void setLastWrite(int i);
 
    dsfmt_t RandomState;
 
 
private:
    int WRITE_FLAG; // Write Flag
    int FLAG; // Flag
    LogPosteriorBase* L_; // Pointer to the LogPosterior base passed in through contstructor
    struct outputInfo{
        string outtype;
        string filename;
        int freq_file;
        int freq_screen;
    } outputinfo_;
    int chainDim_; // Chain dimensions
    double (*logPosterior_)(Array1D<double>&, void *) = NULL; // Pointer to log-posterior function
    void (*fcnAccept_)(void *) = NULL; // Pointer to accept function
    void (*fcnReject_)(void *) = NULL; // Pointer to reject function
    void *postInfo_ = NULL; // Void pointer to the posterior info (data)
    Array2D<double> chcov; // Chain proposal distributions (before the adaptivity starts)
    int seed_; // Random seed for mcmc
 
    virtual double probOldNew(Array1D<double>& a, Array1D<double>& b){return 0.0;}; // Evaluate old|new probabilities and new|old probabilities
    //double evallogMVN_diag(Array1D<double>& x,Array1D<double>& mu,Array1D<double>& sig2); // Evaluate MVN
 
    chainstate currState_; // The current chain state
    chainstate modeState_; // The current MAP state
    Array1D<chainstate> fullChain_; // Array of chain states
 
    //void updateMode(); // Function to update the chain mode
 
    int lastwrite_; // Indicates up to which state
    bool namesPrepend = false;
 
    bool newMode_ = false; // Flag to indicate whether a new mode is found during the last call of runChain, initalized as false
 
    double accRatio_ = -1.0; // Acceptance ratio of the chain, initialized as -1.0
 
    // Flags to indicate whether the corresponding parameters are initialized or not
    bool chaindimInit_ = false;
    bool propcovInit_ = false;
    bool methodInit_ = false;
    bool outputInit_ = false;
 
    bool fcnAcceptFlag_ = false; // Flag that indicates whether the accept function is given or not
    bool fcnRejectFlag_ = false; // Flag that indicates whether the reject function is given or not
 
    Array1D<double> Lower_;
    Array1D<double> Upper_;
 
    Array1D<int> lower_flag_;
    Array1D<int> upper_flag_;
 
};
 
#endif /* UQTKMCMC_H_SEEN */