As I always say, we take pull requests.
These things always sound easy when you just drop one liners into the group. When you get down and wrestle with the guts of what's going on, you'll see things are much more complex. Let's just take the first step at trying to flesh out "take a simple functor".
When you say "simple functor", do you mean one with a templated method like this:
T operator()(const Eigen::Matrix<T, -1, 1>& theta) const;
That is, we assume the functor provides something we can run with autodiff? Or do you want to supply a class instance with two operations, one to do the log density and one to do gradients? Or maybe one method to compute both and return
pair<double, VectorXd> operator()(const Eigen::VectorXd& theta) const;
where the first return element is the value and the second is the gradient? Or maybe for efficiency we pass the gradient in?
double operator()(const Eigen::VectorXd& theta, Eigen::VectorXd& grad) const;
where it returns the log density and sets the gradient?
What do we do about saving intermediate output? If it's a simple functor, where do print statements go in the model that generated the functor? How do things like rejections in the body of the functor get handled and reported? Presumably that would mean wrapping up the model class into a simple functor where all this stuff was plumbed in to the evaluation. What do we then do about saving output such as warnings, rejections, intermediate state, etc.?