base/baseHamiltonian.h

#ifndef _BETA_BASEHAMILTONIAN_

#include <Eigen/Core>
#include <RandomLib/Random.hpp>

// Disable Eigen runtime checks
#define NDEBUG

using namespace Eigen;
using RandomLib::Random;




class baseHamiltonian
{
  
    public:
    
        explicit baseHamiltonian(int dim);
        virtual ~baseHamiltonian() {};
    
        //                   Accessors                  //
    
        double dim() { return mDim; }              

        virtual double T() = 0;                    
        virtual double V() = 0;                    
        virtual const VectorXd& gradV() = 0;       
    
        double H() { return T() + V(); }           
    
        VectorXd& q() { return mQ; }               
        double q(int i) { return mQ(i); }          
        VectorXd& p() { return mP; }               
        double p(int i) { return mP(i); }          
    
        VectorXd& acceptQ() { return mAcceptQ; }   
        VectorXd& rejectQ() { return mRejectQ; }   
    
        virtual void evolveQ(const double epsilon) = 0;  
    
        //                   Mutators                   //

        virtual void prepareEvolution() {};
    
        virtual void beginEvolveP(const double epsilon) = 0;
        
        virtual void finishEvolveP(const double epsilon) = 0;  
        
        virtual void bounceP(const VectorXd& normal) = 0;
    
        virtual void sampleP(Random& random) = 0;
    
        void saveCurrentPoint() { mStoreQ = mQ, mStoreP = mP; }
    
        void restoreStoredPoint() { mQ = mStoreQ, mP = mStoreP; }
        
        void saveAsRejectSample() { mRejectQ = mQ; }
    
        void saveAsAcceptSample() { mAcceptQ = mQ; }
        
        void sampleWindows(const bool accept) { accept ? mQ = mAcceptQ : mQ = mRejectQ; }
    
        void setAverageDecay(const double alpha) { mMovingAlpha = alpha; }
    
        void updateMetroStats(const bool b, const double a);
    
        void clearHistory();
    
        //              Auxiliary Functions             //
    
        bool isNaN();
    
        virtual bool supportViolated() { return false; }
        
        virtual VectorXd& supportNormal() { return mN; }
    
        virtual void checkEvolution(const double epsilon = 1e-6);
    
        double acceptRate() { return mNumAccept / (mNumAccept + mNumReject); }
        
        double movingAcceptRate() { return mAcceptRateBar; }
    
        virtual void displayState();
    
    protected:
    
        const int mDim;       
    
        VectorXd mQ;          
        VectorXd mStoreQ;     
        VectorXd mRejectQ;    
        VectorXd mAcceptQ;    
    
        VectorXd mP;          
        VectorXd mStoreP;     
    
        VectorXd mN;          
    
        VectorXd mGradV;      
    
        // Metropolis parameters
        double mNumAccept;        
        double mNumReject;        
        
        double mAcceptRateBar;    
        double mEffN;             
        double mMovingAlpha;      
    
};

#define _BETA_BASEHAMILTONIAN_
#endif