One redundant gradient evaluation on every iteration

Summary:

The number of log prob gradient evaluations per sample is one greater than the reported n_leapfrog for that sample.
This does not need to be so; the gradient for the starting point was calculated in the previous iteration and has been saved by the sampler.

Description:

Every transition in the sampler begins with initializing this->z_ and sampling the momentum
https://github.com/stan-dev/stan/blob/274a93c2d670edea4918bb80b1f278f55a9de654/src/stan/mcmc/hmc/nuts/base_nuts.hpp#L82-L85
and ends with copying the selected sample to this->z_
https://github.com/stan-dev/stan/blob/274a93c2d670edea4918bb80b1f278f55a9de654/src/stan/mcmc/hmc/nuts/base_nuts.hpp#L201
But stan::services::util::generate_transitions() (the only functions that calls sampler.transition()) just passes around the sample it got from sampler.transition(). What's the point of copying z_ back and forth and re-initializing the gradient when this->z_ already has the correct gradient from the previous iteration?

(Incidentally: why does it call hamiltonian_.sample_p() before calling hamiltonian_.init()? In principle a new point does not have its metric set until init() and sampling a momentum needs the metric. This doesn't really matter for Euclidean metrics, which are constant, but I'd expect it to break the softabs metric. I couldn't test that theory because I don't know how to expose a sampler with softabs in the services.)

Simply deleting the redundant this->hamiltonian_.init(this->z_, logger); causes no change in output and, for example, a 3d gaussian model runs a couple percent faster. That difference is barely distinguishable from noise and would be expected to be even less for a more complex model. As far as performance is concerned this doesn't matter much.
I guess the current design was chosen because it leads to a (marginally simpler?) API that treats the sampler as stateless. I disagree with that; Markov chains are naturally stateful and the sampler API should not hide it.
For example, split the transition() method into three phases.

// new stateful API
virtual void set_sample(sample&, logger&) = 0;
virtual void update_state(logger&) = 0;
virtual sample get_sample() = 0;

// old, stateless API implemented with the stateful API
sample transition(sample& s, logger& logger) {
  this->set_sample(s, logger);
  this->update_state(logger);
  return this->get_sample();
}

Then generate_transitions() does not need to call set_sample().

Current Version:

v2.28.1

Read the referenced sections of src/stan/mcmc/hmc/nuts/base_nuts.hpp and trace stan::services::util::generate_transitions() through sampler.transition(). Check the proposed stateful transition phases and the ordering of hamiltonian_.sample_p() and hamiltonian_.init(). Done means removing the redundant gradient work without changing output, validated with the mentioned 3d Gaussian model.