base/chainBundle.hpp

#ifndef _BETA_CHAINBUNDLE_

#include <vector>
#include <RandomLib/Random.hpp>

using std::vector;
using RandomLib::Random;
                                           

template <typename T>
class chainBundle
{
    
    public:
        
        explicit chainBundle(Random &random);
        ~chainBundle();
        
        //                   Accessors                  //
        
        int nChains() { return (int)mChains.size(); }                  
        T* chain(const int i) { return mChains.at(i); }                
    
        double stepSize() { return mStepSize; }                        
        double nLeapfrog() { return mNumLeapfrog; }                    
    
        double ESS(int i) { return mESS(i); }                          
        double minESS() { return mESS.minCoeff(); }                    
        
        vector<VectorXd>& samples() { return mSamples; }               
    
        //                   Mutators                   //
        
        void setVerbosity(const unsigned int v) { mVerbosity = v; }
        
        void setNumSubsample(const int n) { mNumSubsample = n; }        
        void setNumLeapfrog(const int n) { mNumLeapfrog = n; }          
        void setWindowSize(const int n);                                //   Requires warning, see cpp
        
        void setStepSize(const double s) { mStepSize = s; }             
        void setStepSizeJitter(const double j);                         //   Requires warning, see cpp
        void setStepSizeJerk(const double j);                           //   Requires warning, see cpp
        void setProbStepSizeJerk(const double p);                       //   Requires warning, see cpp

        void setTemperAlpha(const double a) { mTemperAlpha = a; }       
        
        void setAdaptIntTime(double t) { mAdaptIntTime = t; }           
        void setAdaptTargetAccept(double a) { mAdaptTargetAccept = a; } 
        void setAdaptMaxLeapfrog(int n) { mAdaptMaxLeapfrog = n; }      
        
        void storeSamples(const bool store = true) { mStoreSamples = store; }  
        
        void setMaxLagDisplay(int n) { mMaxLagDisplay = n; } 
        void setMaxLagCalc(int n) { mMaxLagCalc = n; }       
                              
        void useExternalMoments(VectorXd& mean, VectorXd& var) { mUseExternalMoments = true; mExternalMean = &mean; mExternalVar = &var; }
    
        void addChain(T* object) { mChains.push_back(object); }
        
        //              Auxiliary Functions             //
        
        void evolve(T* chain, const double epsilon);
        void verboseEvolve(T* chain, const double epsilon);

        void engageAdaptation() { mAdaptStepSize = true; }     
        void disengageAdaptation() { mAdaptStepSize = false; } 
        
        void initAdaptation();
        void updateAdaptation(double metropolisProb);
    
        void transition();
        void transition(int i);
        bool transition(T* chain);
    
        void saveTrajectory(std::ofstream& outputFile, const int i = 0);
    
        void checkIntError(const int i, const int N, const int f = 1);
        
        void seed(const double min, const double max);
        void seed(const int i, const double min, const double max);
        void burn(const int nBurn = 100, const int nCheck = 100, const double minR = 1.1, const bool verbose = true);
        
        void clearSamples();
        void computeSummaryStats(const int b, const bool verbose = true, std::ofstream* outputFile = NULL);
        
    protected:
        
        // Flags
        bool mStoreSamples;       
        bool mAdaptStepSize;      
        bool mRefStats;           
        bool mChar;               
        bool mConstraint;         
        bool mUseExternalMoments; 
        
        // Switches
        unsigned int mVerbosity;  
        
        // Random number generator
        Random mRandom;           
        
        // Hamiltonian Monte Carlo parameters
        int mNumSubsample;        
        int mNumLeapfrog;         
        int mWindowSize;          
        
        double mStepSize;         
        double mStepSizeJitter;   
        double mStepSizeJerk;     
        double mProbStepSizeJerk; 
        
        int mNumFixedPoint;       
        
        double mTemperAlpha;      
    
        double mAdaptTargetAccept; 
        double mAdaptIntTime;      
        int mAdaptMaxLeapfrog;     
        double mAdaptMu;           
        double mAdaptLogEpsBar;    
        double mAdaptHBar;         
        double mAdaptCounter;      
        double mAdaptGamma;        
        double mAdaptKappa;        
        double mAdaptT0;           
        
        // Sample parameters
        int mNumSamples;             
        vector<VectorXd> mSamples;   
        vector<double> mSampleE;     

        int mMaxLagDisplay;          
        int mMaxLagCalc;             
    
        VectorXd mESS;           
    
        VectorXd* mExternalMean;     
        VectorXd* mExternalVar;      

        // Chain containers
        vector<T*> mChains;       
        typename vector<T*>::iterator mIt; 
        
        // Private functions
    
        virtual bool fConstraint(T* chain) { return !(chain->supportViolated()); }
        
        virtual const VectorXd& fNormal(T* chain) { return chain->supportNormal(); }
     
        double fLogSumExp(double a, double b);
    
};

// C++ standard libraries
#include "math.h"
#include <iostream>
#include <iomanip>
#include <fstream>

using std::cout;
using std::endl;
using std::flush;

template <typename T>
chainBundle<T>::chainBundle(Random& random): 
mRandom(random),
mESS(VectorXd::Zero(1))
{
    
    mStoreSamples = true;
    mAdaptStepSize = false;
    mRefStats = false;
    mChar = false;
    mConstraint = true;
    mUseExternalMoments = false;
    
    mVerbosity = 0;
    
    mNumSubsample = 1;
    mNumLeapfrog = 150;
    mWindowSize = 10;
    
    mStepSize = 0.1;
    mStepSizeJitter = 0.05;
    mStepSizeJerk = 0.1;
    mProbStepSizeJerk = 0.0;
    
    mTemperAlpha = 0.0;
    
    mAdaptTargetAccept = 0.65;
    mAdaptIntTime = 10;
    mAdaptMaxLeapfrog = 1000;
    mAdaptMu = 0;
    mAdaptLogEpsBar = 0;
    mAdaptHBar = 0;
    mAdaptCounter = 0;
    mAdaptGamma = 0.05;
    mAdaptKappa = 0.75;
    mAdaptT0 = 10;   
    
    mNumSamples = 0;
    
    mMaxLagDisplay = 10;
    mMaxLagCalc = 50;

    mExternalMean = 0;
    mExternalVar = 0;
    
}



template <typename T>
chainBundle<T>::~chainBundle()
{
    for(mIt = mChains.begin(); mIt != mChains.end(); ++mIt) delete *mIt;
    mChains.clear();
    
    clearSamples();
    
}


template <typename T>
void chainBundle<T>::setWindowSize(const int n)
{
    
    if(n > 0.5 * mNumLeapfrog)
    {
        cout << "chainBundle::setWindowSize() - Reject/Accept window must not be larger"
        << " than half of the total leapfrog steps!" << endl;
        return;
    }
    
    if(n < 2)
    {
        //cout << "chainBundle::setWindowSize() - Reject/Accept window should be larger than two"
        //     << " to avoid proposed windows in violation of the defined constraints!" << endl; 
        //return;
    }
    
    mWindowSize = n;
    
}


template <typename T>
void chainBundle<T>::setStepSizeJitter(const double j)
{
    
    if(j < 0 || j > 1)
    {
        cout << "chainBundle::setStepSizeJitter() - Step size jitter must lie between 0 and 1!" << endl;
        return;
    }
    
    mStepSizeJitter = j;
    
}


template <typename T>
void chainBundle<T>::setStepSizeJerk(const double j)
{
    
    if(j < 0)
    {
        cout << "chainBundle::setStepSizeJerk() - Step size jerk cannot be negative!" << endl;
        return;
    }
    
    mStepSizeJitter = j;
    
}


template <typename T>
void chainBundle<T>::setProbStepSizeJerk(const double p)
{
    
    if(p < 0 || p > 1)
    {
        cout << "chainBundle::setProbStepSizeJerk() - Step size jerk probability must lie between 0 and 1!" << endl;
        return;
    }
    
    mProbStepSizeJerk = p;
    
}


template <typename T>
void chainBundle<T>::clearSamples()
{
    mSamples.clear();
    mSampleE.clear();
    for(mIt = mChains.begin(); mIt != mChains.end(); ++mIt) (*mIt)->clearHistory();
    
}


template <typename T>
void chainBundle<T>::computeSummaryStats(const int b, const bool verbose, std::ofstream* outputFile)
{
    
    cout.precision(6);
    int acceptWidth = 14;
    
    if(mSamples.size() == 0)
    {
        cout << "Cannot compute summary stats without any samples!" << endl;
        return;
    }
    
    if(verbose)
    {
        
        cout << "Displaying summary statistics..." << endl << endl;
        
        cout << "    " << std::setw(5 * acceptWidth / 2) << std::setfill('-') << "" << std::setfill(' ') << endl;
        cout << "    "
             << std::setw(acceptWidth / 2) << std::left << "Chain" 
             << std::setw(acceptWidth) << std::left << "Global"
             << std::setw(acceptWidth) << std::left << "Local"
             << std::endl;
        cout << "    "
             << std::setw(acceptWidth / 2) << std::left << "" 
             << std::setw(acceptWidth) << std::left << "Accept Rate"
             << std::setw(acceptWidth) << std::left << "Accept Rate"
             << std::endl;
        cout << "    " << std::setw(5 * acceptWidth / 2) << std::setfill('-') << "" << std::setfill(' ') << endl;
        
        for(int i = 0; i < mChains.size(); ++i)
        {
            cout << "    "
                 << std::setw(acceptWidth / 2) << std::left << i 
                 << std::setw(acceptWidth) << std::left << (mChains.at(i))->acceptRate()
                 << std::setw(acceptWidth) << std::left << (mChains.at(i))->movingAcceptRate()
                 << endl;
        }
        cout << "    " << std::setw(5 * acceptWidth / 2) << std::setfill('-') << "" << std::setfill(' ') << endl;
        cout << endl;
        
    }
    
    int nSamples = (int)mSamples.size();
    int nBatch = nSamples / b;
    int nLag = nBatch / 2;
    
    nLag = nLag > mMaxLagCalc ? mMaxLagDisplay : nLag;
    int nLagDisplay = nLag > mMaxLagDisplay ? mMaxLagDisplay : nLag;
    
    if(verbose)
    {
        cout << "\tUsing " << nBatch << " batches of " << b << " samples"
        << " (" << nSamples % b << " samples clipped)" << endl << endl;
    }
    
    // Loop over samples
    double batchPotential;
    double sumPotential = 0;
    double sumPotential2 = 0;
    
    long dim = chain(0)->dim();
    if(outputFile) *outputFile << "DIM " << dim << endl << endl;
    
    VectorXd batchSample(dim);
    VectorXd sum(VectorXd::Zero(dim));
    VectorXd sum2(VectorXd::Zero(dim));
    
    MatrixXd autoCorr = MatrixXd::Zero(dim, nLag);
    mESS.resize(dim);
    
    vector<VectorXd> batchSamples;
    
    vector<double>::const_iterator eIt = mSampleE.begin();
    vector<VectorXd>::const_iterator it = mSamples.begin();
    
    // Compute batch expectations
    if(outputFile) *outputFile << "BEGINSAMPLES" << endl;
    
    for(int n = 0; n < nBatch; ++n)
    {
        
        batchPotential = 0;
        batchSample.setZero();
        
        // Loop over intra-batch samples
        for(int m = 0; m < b; ++m, ++eIt, ++it)
        {
            
            // Increment sums
            batchPotential += *eIt;
            batchSample += (*it);
            
            // Write to output file, if desired
            if(outputFile)
            {
                *outputFile << n * b + m << "\t" << *eIt << flush;
                for(int i = 0; i < dim; ++i)
                {
                    *outputFile << "\t" << (*it)(i) << flush;
                }
                *outputFile << endl;
            }
        }
        
        // Normalize
        batchPotential /= b;
        sumPotential += batchPotential;
        sumPotential2 += batchPotential * batchPotential;
        
        batchSample /= (double)b;
        batchSamples.push_back(VectorXd(batchSample));
        
        sum += batchSample;
        sum2 += batchSample.array().square().matrix();
        
    }
    
    if(outputFile) *outputFile << "ENDSAMPLES" << endl << endl;
    
    // Normalize batch expectations
    sumPotential /= (double)nBatch;
    sumPotential2 = sumPotential2 / (double)nBatch - sumPotential * sumPotential;
    
    sum /= (double)nBatch;
    
    sum2 /= (double)nBatch; 
    sum2 -= sum.array().square().matrix();
    sum2 *= b / nSamples; // Correct asymptotic batch variance to asymptotic chain variance
    
    // Compute autocorrelations
    VectorXd gamma(dim);
    VectorXd one(dim);
    VectorXd two(dim);
    
    if(!mUseExternalMoments) mExternalMean = &sum;
    
    for(int k = 0; k < nLag; ++k)
    {
        
        gamma = autoCorr.col(k);
        
        vector<VectorXd>::iterator itOne = batchSamples.begin();
        
        vector<VectorXd>::iterator itTwo = itOne;
        for(int j = 0; j < k; ++j) ++itTwo;
        
        for(int n = 0; n < nBatch - k; ++n, ++itOne, ++itTwo)
        {
            one = *itOne - *mExternalMean;
            two = *itTwo - *mExternalMean;
            gamma += one.cwiseProduct(two);
        }
        
        gamma /= (double)nBatch;
        
        // Compute autocorrelation
        if(k == 0) 
        {
            for(int i = 0; i < dim; ++i) autoCorr(i, k) = gamma(i);
            continue;
        }
        
        one = gamma.cwiseQuotient(autoCorr.col(0));
        
        if(mUseExternalMoments) one = gamma.cwiseQuotient(*mExternalVar);
        else                    one = gamma.cwiseQuotient(autoCorr.col(0));
     
        for(int i = 0; i < dim; ++i) autoCorr(i, k) = one(i);
        
    }
    
    autoCorr.col(0).setOnes();
    
    // Compute effective sample size for each variable
    for(int i = 0; i < dim; ++i)
    {
        
        mESS(i) = 0;
        
        double gamma = 0;
        double Gamma = 0;
        double oldGamma = 1;
        
        for(int k = 1; k < nLag; ++k)
        {
            
            gamma = autoCorr(i, k);
            
            // Even k, finish computation of Gamma
            if(k % 2 == 0)
            {
                
                Gamma += gamma;
                
                // Convert to initial monotone sequence (IMS)
                if(Gamma > oldGamma) Gamma = oldGamma;
                
                oldGamma = Gamma;
                
                // Terminate if necessary
                if(Gamma <= 0) break;
                
                // Increment autocorrelation sum
                mESS(i) += Gamma;
                
            }
            // Odd k, begin computation of Gamma
            else
            {
                Gamma = gamma;   
            }
            
            //ess(i) += (1.0 - (double)k / (double)nBatch) * (*(autoCorr[k]))[i];
            
        }
        
        mESS(i) = (double)nBatch / (1.0 + 2 * mESS(i) );
        
    }
    
    // Display results
    if(verbose)
    {
        cout << "\tV ~ " << sumPotential << " +/- " << sqrt(sumPotential2) << endl;
        cout << "\tSmallest Effective Sample Size is " << minESS() << endl << endl;
    }
    
    double whiteNoise = 2.0 / sqrt((double)nBatch);
    
    if(outputFile) 
    {
        
        *outputFile << "WHITENOISE " << whiteNoise << endl << endl;
        *outputFile << "BEGINAUTOCORR" << endl;
        
        for(int k = 0; k < nLag; ++k)
        {
            
            *outputFile << k << flush;
            
            for(int i = 0; i < dim; ++i)
            {
                *outputFile << "\t" << autoCorr(i, k) << flush;
            }
            
            *outputFile << endl;
            
        }
        
        *outputFile << "ENDAUTOCORR" << endl;
        
    }
    
    if(verbose)
    {
        
        for(int i = 0; i < dim; ++i)
        {
            
            cout << "\t" << std::setw(4 * acceptWidth) << std::setfill('-') << "" << std::setfill(' ') << endl;
            cout << "\tVariable " << i << endl;
            cout << "\t" << std::setw(4 * acceptWidth) << std::setfill('-') << "" << std::setfill(' ') << endl;
            cout << endl;
            
            cout << "\t\tMean = " << sum[i] << endl;
            cout << "\t\tMonte Carlo Standard Error = " << sqrt(sum2[i]) << endl;
            cout << "\t\tEffective Sample Size = " << mESS(i) << endl << endl;
            
            cout << endl << "\t\tAutocorrelation:" << endl;
            cout << "\t\t\tWhite noise correlation = " << whiteNoise << endl << endl;
            
            
            for(int k = 0; k < nLagDisplay; ++k)
            {
                double g = autoCorr(i, k);
                cout << "\t\t\tgamma_{" << k << "} = " << g << flush;
                if(g > whiteNoise && k > 0) cout << " (Larger than expected fluctuation from white noise!)" << flush;
                cout << endl;
            }
            
            cout << endl;
            
        }
        
        cout << "\t" << std::setw(4 * acceptWidth) << std::setfill('-') << "" << std::setfill(' ') << endl;
        
    }
    
    // Clean up
    batchSamples.clear();
    
    return;
    
}


template <typename T>
void chainBundle<T>::evolve(T* chain, const double epsilon)
{
    
    // Initial half momentum evolution
    if(mConstraint) chain->beginEvolveP(0.5 * epsilon);
    
    // Full spatial evolution
    chain->evolveQ(epsilon);
    
    // Check constraint
    mConstraint = fConstraint(chain);
    
    // Final half momentum evolution
    if(mConstraint)  chain->finishEvolveP(0.5 * epsilon);
    else             chain->bounceP(fNormal(chain));
    
}


template <typename T>
void chainBundle<T>::verboseEvolve(T* chain, const double epsilon)
{
    
    //cout << "----------------------------------------------------" << endl;
    
    std::cout.precision(6);
    int width = 14;
    int nColumn = 4;
    
    std::cout << "Verbose Hamiltonian Evolution, Step Size = " << epsilon << ":" << std::endl;
    std::cout << "    " << std::setw(nColumn * width) << std::setfill('-') << "" << std::setfill(' ') << std::endl;
    std::cout << "    "
    << std::setw(width) << std::left << "Poisson"
    << std::setw(width) << std::left << "Initial" 
    << std::setw(width) << std::left << "Current"
    << std::setw(width) << std::left << "DeltaH"
    << std::endl;
    std::cout << "    "
    << std::setw(width) << std::left << "Operator"
    << std::setw(width) << std::left << "Hamiltonian" 
    << std::setw(width) << std::left << "Hamiltonian"
    << std::setw(width) << std::left << "/ Stepsize^{2}"
    << std::endl;
    std::cout << "    " << std::setw(nColumn * width) << std::setfill('-') << "" << std::setfill(' ') << std::endl;
    
    double H0 = chain->H();
    
    // Initial half momentum evolution
    if(mConstraint) chain->beginEvolveP(0.5 * epsilon);
    
    double H1 = chain->H();
    
    std::cout << "    "
    << std::setw(width) << std::left << "hat{V}/2"
    << std::setw(width) << std::left << H0
    << std::setw(width) << std::left << H1
    << std::setw(width) << std::left << (H1 - H0) / (epsilon * epsilon)
    << std::endl;
    
    // Full spatial evolution
    chain->evolveQ(epsilon);
    
    double H2 = chain->H();
    
    std::cout << "    "
    << std::setw(width) << std::left << "hat{T}"
    << std::setw(width) << std::left << H0
    << std::setw(width) << std::left << H2
    << std::setw(width) << std::left << (H2 - H0) / (epsilon * epsilon)
    << std::endl;
    
    // Check constraint
    mConstraint = fConstraint(chain);
    
    // Final half momentum evolution
    if(mConstraint)  chain->finishEvolveP(0.5 * epsilon);
    else             chain->bounceP(fNormal(chain));
    
    double H3 = chain->H();
    
    std::cout << "    "
    << std::setw(width) << std::left << "hat{V}/2"
    << std::setw(width) << std::left << H0
    << std::setw(width) << std::left << H3
    << std::setw(width) << std::left << (H3 - H0) / (epsilon * epsilon)
    << std::endl;
    
    std::cout << "    " << std::setw(nColumn * width) << std::setfill('-') << "" << std::setfill(' ') << std::endl;
    
}


template <typename T>
void chainBundle<T>::initAdaptation()
{
    
    mAdaptMu = log(mStepSize / (1 - mStepSize) );
    mAdaptLogEpsBar = log(mStepSize / (1 - mStepSize) );
    
    mAdaptHBar = 0;
    mAdaptCounter = 0;
    
    mNumLeapfrog = (int)(mAdaptIntTime / mStepSize);
    mNumLeapfrog = mNumLeapfrog < 1 ? 1 : mNumLeapfrog;
    mNumLeapfrog = mNumLeapfrog > mAdaptMaxLeapfrog ? mAdaptMaxLeapfrog : mNumLeapfrog;
    
}


template <typename T>
void chainBundle<T>::updateAdaptation(double metropolisProb)
{
    
    metropolisProb = metropolisProb > 1 ? 1 : metropolisProb;
    
    // Update averaging parameters
    ++mAdaptCounter;
    
    double delta = 1.0 / (mAdaptCounter + mAdaptT0);
    
    mAdaptHBar = (1.0 - delta) * mAdaptHBar + delta * (mAdaptTargetAccept - metropolisProb);
    
    double logEps = mAdaptMu - mAdaptHBar * sqrt(mAdaptCounter) / mAdaptGamma;
    delta = exp( - mAdaptKappa * log( mAdaptCounter ) );
    
    mAdaptLogEpsBar = (1.0 - delta) * mAdaptLogEpsBar + delta * logEps;
    
    // Then the stepsize
    mStepSize = 1.0 / ( 1.0 + exp(-mAdaptLogEpsBar) );
    
    // And the number of leapfrog steps to maintain a constant integration time
    mNumLeapfrog = (int)(mAdaptIntTime / mStepSize);
    mNumLeapfrog = mNumLeapfrog < 1 ? 1 : mNumLeapfrog;
    mNumLeapfrog = mNumLeapfrog > mAdaptMaxLeapfrog ? mAdaptMaxLeapfrog : mNumLeapfrog;
    
}


template <typename T>
void chainBundle<T>::transition()
{
    for(mIt = mChains.begin(); mIt != mChains.end(); ++mIt) transition(*mIt);
}


template <typename T>
void chainBundle<T>::transition(int i)
{
    
    if(i < 0 || i >= mChains.size())
    {
        cout << "chainBundle::transition() - Bad chain index " << i << "!" << endl;
        return;
    }
    
    transition(mChains.at(i));
    
    return;
    
}


template <typename T>
void chainBundle<T>::saveTrajectory(std::ofstream& outputFile, const int i)
{
    
    outputFile << mStepSize << "\t" << mNumLeapfrog << endl << endl;
    
    chain(i)->prepareEvolution();
    chain(i)->sampleP(mRandom);
    
    outputFile << 0 << "\t" << chain(i)->H() << "\t" << chain(i)->T() << "\t" << chain(i)->V() << flush;
    
    for(int j = 0; j < chain(i)->dim(); ++j) outputFile << "\t" << chain(i)->q()(j) << flush;
    for(int j = 0; j < chain(i)->dim(); ++j) outputFile << "\t" << chain(i)->p()(j) << flush;  
    outputFile << endl;
    
    for(int n = 0; n < mNumLeapfrog; ++n) 
    {
        
        evolve(chain(i), mStepSize);
        
        outputFile << n + 1 << "\t" << chain(i)->H() << "\t" << chain(i)->T() << "\t" << chain(i)->V() << flush;
        
        for(int j = 0; j < chain(i)->dim(); ++j) outputFile << "\t" << chain(i)->q()(j) << flush;
        for(int j = 0; j < chain(i)->dim(); ++j) outputFile << "\t" << chain(i)->p()(j) << flush;  
        outputFile << endl;
        
    }
    
    outputFile.close();
    
    return;
    
}


template <typename T>
bool chainBundle<T>::transition(T* chain)
{
    
    // Check that the chain was properly seeded
    if(mNumSamples == 0)
    {
        
        mConstraint = fConstraint(chain);
        
        if(!mConstraint)
        {
            cout << "chainBundle::transition() - Initial chain parameters in violation "
            << "of constraints, aborting the transition!" << endl;
            return false;
        }
        
    }
    
    bool accept = true;
    
    for(int n = 0; n < mNumSubsample; ++n)
    {
        
        // Add random jitter to the step size to avoid closed loops, and
        // an occassional jerk to allow the chain to move through narrow valleys
        double stepSize = mStepSize;
        
        if(mRandom.Prob<double>(mProbStepSizeJerk)) 
        {
            stepSize *= mStepSizeJerk;
        }
        else
        {
            stepSize *= 1.0 + mStepSizeJitter * (2.0 * mRandom.FloatU() - 1.0);
        }
        
        // Prepare the chain
        chain->prepareEvolution();
        
        // Sample momenta conditioned on the current position of the chain
        chain->sampleP(mRandom);
        
        //        Sample from the reject window         //
        
        // Immediately accept initial point as current sample
        chain->saveAsRejectSample();
        double logSumRejectProb = -chain->H();
        
        if(mWindowSize > 0)
        {
            
            // Backwards evolution
            stepSize *= -1;
            
            chain->saveCurrentPoint();
            int s = mRandom.Integer(mWindowSize);
            
            for(int t = 0; t < s; ++t)
            {
                
                evolve(chain, stepSize);
                
                // Accept or reject updating the sample, with constraint violations immediately rejected 
                if(mConstraint) 
                {
                    double Hplus = chain->H();
                    logSumRejectProb = fLogSumExp(logSumRejectProb, -Hplus);
                    if( mRandom.Prob<double>( exp(-Hplus - logSumRejectProb) ) ) chain->saveAsRejectSample();
                }
                
            }
            
            // Forwards evolution
            stepSize *= -1;
            
            chain->restoreStoredPoint();
            
            for(int t = 0; t < mWindowSize - s - 1; ++t)
            {
                
                evolve(chain, stepSize);
                
                // Accept or reject updating the sample, with constraint violations immediately rejected 
                if(mConstraint) 
                {
                    double Hplus = chain->H();
                    logSumRejectProb = fLogSumExp(logSumRejectProb, -Hplus);
                    if( mRandom.Prob<double>( exp(-Hplus - logSumRejectProb) ) ) chain->saveAsRejectSample();
                }
                
            }
            
        }
        
        //       Evolve for mNumLeapfrog iterations     //
        
        if(mVerbosity >= 2)
        {
            cout << "Beginning leapfrog:" << endl;
            cout << endl;
            chain->displayState();
            cout << endl;
        }
        
        if(mTemperAlpha) chain->p() *= sqrt(mTemperAlpha);
        
        
        double m = 0.5 * (mNumLeapfrog - 2 * mWindowSize);
        
        for(unsigned int t = 0; t < mNumLeapfrog - 2 * mWindowSize; ++t)
        {
            
            // Immediately reject if a NaN shows up
            if(chain->isNaN()) 
            {
                if(mVerbosity >= 1) cout << "Immediately rejecting a NaN state..." << endl;
                break;
            }
            
            // Evolve one step
            evolve(chain, stepSize);
            
            // Temper the momenta
            if(mTemperAlpha)
            {
                if(t < m)       chain->p() *= sqrt(mTemperAlpha);
                else if(t >= m) chain->p() /= sqrt(mTemperAlpha);
            }
            
        }
        
        if(mVerbosity >= 2)
        {
            cout << "Ending leapfrog:" << endl;
            cout << endl;
            chain->displayState();
            cout << endl;
        }
        
        //        Sample from the accept window         //
        
        double logSumAcceptProb = 0;
        bool goodAcceptSample = false;
        
        if(!chain->isNaN())
        {
            
            // Immediately accept current point, provided the constraint is not violated            
            if(mConstraint) 
            {
                chain->saveAsAcceptSample();
                logSumAcceptProb = -chain->H();
                goodAcceptSample = true;
            }
            
            if(mWindowSize > 0)
            {
                
                for(int t = 0; t < mWindowSize - 1; ++t)
                {
                    
                    evolve(chain, stepSize);
                    
                    // Accept or reject updating the sample, with constraint violations immediately rejected 
                    if(mConstraint) 
                    {
                        
                        // Immediately reject if a NaN shows up
                        if(chain->isNaN()) 
                        {
                            if(mVerbosity >= 1) cout << "Immediately rejecting a NaN state..." << endl;
                            break;
                        }
                        
                        double Hplus = chain->H();
                        
                        // Immediately accept if no sample has yet been accepted
                        if(!goodAcceptSample)
                        {
                            logSumAcceptProb = -Hplus;
                            chain->saveAsAcceptSample();
                            goodAcceptSample = true;
                        }
                        else
                        {
                            logSumAcceptProb = fLogSumExp(logSumAcceptProb, -Hplus);
                            if( mRandom.Prob<double>( exp(-Hplus - logSumRejectProb) ) ) chain->saveAsAcceptSample();
                        }
                        
                    }
                    
                }
                
            }
            else
            {
                
                if(!mConstraint)
                {
                    chain->saveAsAcceptSample();
                    goodAcceptSample = false;
                }
                
            }
            
        }
        
        //   Sample between reject and accept windows   //
        
        double metropolisProb = 0;
        if(goodAcceptSample) metropolisProb = exp(logSumAcceptProb - logSumRejectProb);
        if(metropolisProb != metropolisProb) metropolisProb = 0;
        
        if(metropolisProb > 1) accept = true;
        else                   accept = mRandom.Prob<double>(metropolisProb);
        
        chain->sampleWindows(accept);
        
        if(mVerbosity >= 2)
        {
            if(accept) cout << "Accepting proposed sample with" << endl;
            else       cout << "Rejecting proposed sample with" << endl;
            
            cout << "    log(p_{accept window}) = " << logSumAcceptProb << endl;
            cout << "    log(p_{reject window}) = " << logSumRejectProb << endl;
            cout << "    p(accept window) = " << metropolisProb << endl << endl;
            
        }
        
        chain->updateMetroStats(accept, metropolisProb > 1 ? 1.0 : metropolisProb);
        
        if(mAdaptStepSize) updateAdaptation(metropolisProb);
        
    }
    
    // Store sample if desired
    if(mStoreSamples) 
    {
        mSamples.push_back(VectorXd(chain->q()));
        mSampleE.push_back(chain->V());
    }
    
    ++mNumSamples;
    
    return accept;
    
}


template <typename T>
void chainBundle<T>::checkIntError(int i, int N, int f)
{
    
    if(i < 0 || i >= mChains.size())
    {
        cout << "chainBundle::checkIntError() - Bad chain index " << i << "!" << endl;
        return;
    }
    
    // Prepare chain
    mChains.at(i)->prepareEvolution();
    mChains.at(i)->sampleP(mRandom);
    
    mChains.at(i)->displayState();
    
    double initH = mChains.at(i)->H();
    
    // Display Header
    std::cout.precision(6);
    int width = 12;
    int nColumn = 4;
    
    std::cout << "Displaying the integration error of chain " << i << ":" << std::endl;
    std::cout << "    " << std::setw(nColumn * width) << std::setfill('-') << "" << std::setfill(' ') << std::endl;
    std::cout << "    "
    << std::setw(width) << std::left << "Step" 
    << std::setw(width) << std::left << "H_{0}"
    << std::setw(width) << std::left << "H_{i}"
    << std::setw(width) << std::left << "H_{0} - H_{i}"
    << std::endl;
    std::cout << "    " << std::setw(nColumn * width) << std::setfill('-') << "" << std::setfill(' ') << std::endl;
    
    // Evolve and display results
    for(int n = 0; n < N; ++n)
    {
        
        evolve(mChains.at(i), mStepSize);
        
        if(!(n % f))
        {
            std::cout << "    "
            << std::setw(width) << std::left << n
            << std::setw(width) << std::left << initH
            << std::setw(width) << std::left << mChains.at(i)->H()
            << std::setw(width) << std::left << initH -  mChains.at(i)->H()
            << std::endl;
        }
        
    }
    
    std::cout << "    " << std::setw(nColumn * width) << std::setfill('-') << "" << std::setfill(' ') << std::endl;
    std::cout << std::endl;
    
}


template <typename T>
void chainBundle<T>::seed(const double min, const double max)
{
    for(int i = 0; i < mChains.at(0)->dim(); ++i) seed(i, min, max);
}


template <typename T>
void chainBundle<T>::seed(const int i, const double min, const double max)
{
    
    if(i < 0 || i >= mChains.at(0)->dim())
    {
        cout << "chainBundle::seed() - Bad parameter index!" << endl;
        return;
    }
    
    for(mIt = mChains.begin(); mIt != mChains.end(); ++mIt) 
    {
        (*mIt)->q()(i) = (max - min) * mRandom.FloatU<double>() + min;
    }
    
}


template <typename T>
void chainBundle<T>::burn(const int nBurn, const int nCheck, const double minR, const bool verbose)
{
    
    // Convergence storage
    long dim = mChains.at(0)->dim();
    double m = (double)mChains.size();
    double n = (double)nCheck;
    
    VectorXd R(dim);   // Gelman-Rubin statistics
    
    VectorXd* chainSum[(int)m];
    for(int i = 0; i < m; ++i) chainSum[i] = new VectorXd(dim);
    VectorXd ensembleSum(dim);
    VectorXd ensembleSum2(dim);
    VectorXd chainSum2(dim);
    
    bool burning = true;
    
    int nPerLine = 5;
    int nLines = (int)dim / nPerLine;
    
    cout.precision(6);
    int width = 12;
    int totalWidth = width + 8;
    int borderWidth = ((int)dim < nPerLine ? (int)dim : nPerLine) * totalWidth + 4;
    
    int acceptWidth = 14;
    
    //              Burn in Markov Chain            //
    
    // Temporarily store samples
    bool sampleFlag = mStoreSamples;
    mStoreSamples = false;
    
    // Burn through initial samples
    if(verbose) cout << endl << "Burn, baby, burn" << flush;
    
    for(int i = 0; i < nBurn; ++i) 
    {
        transition();
        if(verbose) if(i % 10 == 0) cout << "." << flush;
    }
    
    if(verbose) cout << endl << endl;
    if(verbose) cout << "Computing diagnostics:" << endl;
    if(verbose) cout << "    " << std::setw(borderWidth) << std::setfill('-') << "" << std::setfill(' ') << endl;
    
    while(burning)
    {
        
        ensembleSum.setZero();
        ensembleSum2.setZero();
        for(int i = 0; i < m; ++i) (chainSum[i])->setZero();
        chainSum2.setZero();
        
        // Burn through diagnostic samples
        for(int i = 0; i < nCheck; ++i) 
        {
            
            transition();
            
            mIt = mChains.begin();
            for(int j = 0; j < m; ++j, ++mIt)
            {
                
                VectorXd& q = (*mIt)->q();
                
                ensembleSum += q;
                ensembleSum2 += q.cwiseAbs2();
                
                *(chainSum[j]) += q;
                
            }
            
        }
        
        // Compute Gelman-Rubin statistic for each parameter
        double a = (n - 1) / n;
        double b = m / (m - 1);
        
        for(int j = 0; j < m; ++j)
        {
            chainSum2 += (chainSum[j])->cwiseAbs2();
        }
        
        ensembleSum = ensembleSum.cwiseAbs2();
        ensembleSum = ensembleSum.cwiseQuotient(chainSum2);
        ensembleSum /= m;
        ensembleSum = VectorXd::Ones(dim) - ensembleSum;
        
        ensembleSum2 = ensembleSum2.cwiseQuotient(chainSum2);
        ensembleSum2 *= n;
        ensembleSum2 = VectorXd::Ones(dim) - ensembleSum2;
        
        R = a * (VectorXd::Ones(dim) - b * ensembleSum.cwiseQuotient(ensembleSum2));
        
        burning = false;
        for(int i = 0; i < dim; ++i) burning |= (R[i] > minR);
        
        // Failure output
        if(verbose && burning)
        {
            cout << "    Diagnostic test failed," << endl;
            
            int k = 0;
            for(int i = 0; i < nLines; ++i)
            {
                
                cout << "        " << flush;
                for(int j = 0; j < nPerLine; ++j, ++k)
                {
                    cout << "R_{" << k << "} = " << std::setw(width) << std::left << R[k] << flush;
                }
                
                cout << endl;
                
            }
            
            if(dim % nPerLine > 0)
            {
                
                cout << "        " << flush;
                for(int i = 0; i < dim % nPerLine; ++i, ++k)
                {
                    cout << "R_{" << k << "} = " << std::setw(width) << std::left << R[k] << flush;
                }
                cout << endl;
                
            }
            
            cout << "    " << std::setw(borderWidth) << std::setfill('-') << "" << std::setfill(' ') << endl;
            
        }
        
    }
    
    // Display convergence details if desired
    if(verbose)
    {
        
        cout << "    Markov chains converged with" << endl;
        
        int k = 0;
        for(int i = 0; i < nLines; ++i)
        {
            
            cout << "        " << flush;
            for(int j = 0; j < nPerLine; ++j, ++k)
            {
                cout << "R_{" << k << "} = " << std::setw(width) << std::left << R[k] << flush;
            }
            
            cout << endl;
            
        }
        
        if(dim % nPerLine > 0)
        {
            
            cout << "        " << flush;
            for(int i = 0; i < dim % nPerLine; ++i, ++k)
            {
                cout << "R_{" << k << "} = " << std::setw(width) << std::left << R[k] << flush;
            }
            cout << endl;
            
        }
        
        cout << "    " << std::setw(borderWidth) << std::setfill('-') << "" << std::setfill(' ') << endl;
        
        cout << endl;
        
        cout << "    " << std::setw(5 * acceptWidth / 2) << std::setfill('-') << "" << std::setfill(' ') << endl;
        cout << "    "
        << std::setw(acceptWidth / 2) << std::left << "Chain" 
        << std::setw(acceptWidth) << std::left << "Global"
        << std::setw(acceptWidth) << std::left << "Local"
        << std::endl;
        cout << "    "
        << std::setw(acceptWidth / 2) << std::left << "" 
        << std::setw(acceptWidth) << std::left << "Accept Rate"
        << std::setw(acceptWidth) << std::left << "Accept Rate"
        << std::endl;
        cout << "    " << std::setw(5 * acceptWidth / 2) << std::setfill('-') << "" << std::setfill(' ') << endl;
        
        for(int i = 0; i < mChains.size(); ++i)
        {
            cout << "    "
            << std::setw(acceptWidth / 2) << std::left << i 
            << std::setw(acceptWidth) << std::left << (mChains.at(i))->acceptRate()
            << std::setw(acceptWidth) << std::left << (mChains.at(i))->movingAcceptRate()
            << endl;
        }
        cout << "    " << std::setw(5 * acceptWidth / 2) << std::setfill('-') << "" << std::setfill(' ') << endl;
        cout << endl;
        
    }
    
    mChar = true;
    
    // Return sample flag
    mStoreSamples = sampleFlag;
    
    return;
    
}


template <typename T>
double chainBundle<T>::fLogSumExp(double a, double b)
{
    
    if(a > b) 
    {
        return a + log( 1.0 + exp(b - a) );
    }
    else
    {      
        return b + log( 1.0 + exp(a - b) );
    }
    
}

#define _BETA_CHAINBUNDLE_
#endif