leapfrog.hpp

Go to the documentation of this file.

#pragma once
#include <fastad>
 
namespace ppl {
namespace mcmc {
 
template <class ADExprType
        , class MatType>
inline double reset_autodiff(ADExprType& ad_expr, 
                             Eigen::MatrixBase<MatType>& adjoints,
                             Eigen::MatrixBase<MatType>& tp_adjoints)
{
    adjoints.setZero();
    tp_adjoints.setZero();
    return ad::autodiff(ad_expr);
}
 
template <class ADExprType
        , class MatType
        , class MomentumHandlerType>
inline double leapfrog(ADExprType& ad_expr,
                       Eigen::MatrixBase<MatType>& theta,
                       Eigen::MatrixBase<MatType>& theta_adj,
                       Eigen::MatrixBase<MatType>& tp_adj,
                       Eigen::MatrixBase<MatType>& r,
                       const MomentumHandlerType& m_handler,
                       double epsilon,
                       bool reuse_adj)
{
    if (!reuse_adj) {
        reset_autodiff(ad_expr, theta_adj, tp_adj);
    }
    const double half_step = epsilon/2.;
    r += half_step * theta_adj;
 
    theta += epsilon * m_handler.dkinetic_dr(r);
 
    const double new_potential = 
        -reset_autodiff(ad_expr, theta_adj, tp_adj);
    r += half_step * theta_adj;
 
    return new_potential;
}
 
} // namespace mcmc
} // namespace ppl