finite_diff_gradient_auto.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUNCTOR_FINITE_DIFF_GRADIENT_AUTO_HPP
#define STAN_MATH_PRIM_FUNCTOR_FINITE_DIFF_GRADIENT_AUTO_HPP
 
#include <stan/math/prim/fun/Eigen.hpp>
#include <stan/math/prim/fun/finite_diff_stepsize.hpp>
#include <cmath>
 
namespace stan {
namespace math {
 
template <typename F, typename VectorT, typename GradVectorT,
          typename ScalarT = return_type_t<VectorT>>
void finite_diff_gradient_auto(const F& f, VectorT&& x, ScalarT& fx,
                               GradVectorT& grad_fx) {
  using EigT = Eigen::Matrix<ScalarT, -1, 1>;
  static constexpr int h_scale[6] = {3, 2, 1, -3, -2, -1};
  static constexpr int mults[6] = {1, -9, 45, -1, 9, -45};
 
  fx = f(x);
  grad_fx.resize(x.size());
  Eigen::Map<EigT> grad_map(grad_fx.data(), grad_fx.size());
 
  grad_map = EigT::NullaryExpr(x.size(), [&f, &x](Eigen::Index i) {
    double h = finite_diff_stepsize(value_of_rec(x[i]));
    ScalarT delta_f = 0;
    for (int j = 0; j < 6; ++j) {
      auto x_temp
          = EigT::NullaryExpr(x.size(), [&x, &i, &h, &j](Eigen::Index k) {
              return k == i ? x[i] + h * h_scale[j] : x[k];
            });
      delta_f += f(std::move(x_temp)) * mults[j];
    }
    return delta_f / (60 * h);
  });
}
 
}  // namespace math
}  // namespace stan
#endif