In order to avoid the extensive code duplication that is currently spread across the interfaces we need the core Stan API to start providing analysis routes, including the calculation of effective sample size, rhat, divergences, E-FMI, etc. This will require some reorganization of the stan source.
First and foremost we need to decide on the exact format for I/O, in particular how sample output must be input to the API and how messages and error status is returned by the API. Much of this can be cribbed from the recent discussion about exposing C++ effective sample size functions .
Then we need to reorganize the stan source. Namely we need all analysis code collected together (in particular, effective sample size and rhat functionality should be removed from the chains class as the OO structure is not needed). The combined functionality then needs to be exposed via appropriate services.
Given the partial steps that are being discussed we might as well do this systematically.
I’ll add non-MCMC analysis to this system. Now we have only ADVI, but we may have something else later, but all distributional approximations can use the new Pareto diagnostic. Given draws theta^{(s)} from any distributional approximation q, and log(q(theta^{(s)})) and log(p(theta^{(s)})) we can compute Pareto k (khat) and effective sample size (N_eff) estimates. These should be reported in every interface similarly as Rhat and N_eff are reported for HMC/NUTS.
Yes, once we’ve established conventions for the routes we can add others for other algorithms. We will, however, have to be careful to either enforce the necessary inputs or have well-defined error messages (and return codes) when the inputs are not well-formed for a particular analysis.
I think that’s a good idea. I also think the core implementations could go into stan-dev/math, with only the interfaces and I/O dealt with at the Stan level.
I’m most concerned with getting the C++ APIs in better shape for all these functions. We’ve been discussing the effective sample size calcs in another thread and thought vector<const double*> would be general enough to avoid big copies and accomodate both std::vector array types and Eigen::Matrix vector types.
Currently, I/O isn’t handled at the Stan level—it’s just delegated to the interfaces. RStan has the ability to write out output that CmdStan’s analysis tools can read and also write out data that CmdStan can read. It can also read CmdStan data since that’s just an R format. Are you going to be suggesting restarting the JSON/protobuf discussions?
No, I am suggesting that we define an I/O convention for how the interfaces have to wrap the samples and diagnostic info from a fit object into std::vectors or POD arrays to be analyzed by routes and then what they expect in return (probably strings and error codes if not numerical results like ESS and rhat).
Bonus points if the convention minimized excess copying on the C++ side.
None of this would be exposed to users – the interfaces would call these routes internally.
I am hesitant to put this stuff into math for the same reason I don’t think the Welford stuff really belongs there – the math library is essentially an autodiff library, not a general purpose math library. This stuff would be just as well in a separate auxiliary math library.
That’s a nice way to think about it. Effective sample size estimators and R-hat are more like the sampling algorithms—subject to change as we learn more. So it makes sense to leave those with the algorithms. Plus, I think it makes sense to have you manage them.
If we ever get around to adding derivatives to FFT, that’ll go in the math lib.
Welford’s algorithm is one of those oddball things where it’s at an API level of an accumulator rather than being like our ordinary math functions. So I’m OK with that going either way. I could see some of the algorithms that do need to be in the math library, like calculating sample variance, needing it, though. If that happens, then whatever they use in common should be in the math lib.
Yes, this is what incited my post. We shouldn’t do this piecemeal but put something systematic together on which we can add all of the other diagnostic calculations.
But how would we get the variable/diagnostic x draw x chain shape into a single vector? If we used vectors of vectors then we couldn’t guarantee memory locality for map. We’re probably have to do a copy at least once to put everything together, but let’s figure out the signature that will make it easiest for the interfaces to pack everything together and for the C++ to unpack into Eigen objects.
Chains are stored in separate holders in Python and R. And the variable/diagnostic \times draw is stored using memory-contiguous column-major-ordered holders (e.g., numpy array in Python) by the interfaces so those can be referenced with a double*.
Is there a place where we operate on a Chain x Draw object any other way than as a vector for each chain? I don’t see a need for Variable x Chain x Draw.
If you have vector<double*>, it’s easy to map each entry to an Eigen vector or row vector. So if the answer to the above is that all the underlying operations use Eigen vectors, this is the most economical way to pass things that doesn’t implicate Eigen on the interface side, which is what Allen was requesting.
For the ESS calculations as their are implemented we’d need matrices with shape draw x integration – so can the interfaces guarantee that the info from each chain points to contiguous memory? If so then we’d be looking at a signature like
error_code_type diagnostic_function(std::vector<double*>&,
std::vector<std::string>&,
int n_var, int n_draw, writer& writer);
I’m not sure we need the writers for ess. They add a lot of complexity.
Could we use exceptions instead? Or return a std::pair with the first
value being the result and the second value being a string containing
the error message (or empty if no error). (This is a pattern Go uses.)
Yes, was already thinking that exceptions could handle exceptional events. I don’t think ess has any output other than the return (we don’t need any writers), but other diagnostic functions might?
Yes, the double* for each chain would be contiguous, but they wouldn’t be contiguous with each other.
I’d much prefer exceptions if that’s workable from the interface perspectives.
If there’s an error, the value won’t make sense. Returning both the value and error message is then clunky as only one will ever get used. This error code return is also the tradition in C with functions with signatures of the form int foo(..., result_t* result) where the return is a return code indicating error or not and the function setting result with the actual value calculated.
Unfortunately the conversation has diverged a bit here. Let me try to bring it back a bit.
We do not throw exceptions to the interfaces for the algorithm API. This thread is not about redoing the API but adding additional routes for common diagnostic calculations. If there is interest in changing the client relationship and how error information is propagated then a new thread should be started.
Diagnostic functions return various statuses depending on the thresholds we set (and possibly modifiable by the user). Okay, Warning, and Error messages do not indicate that the function has failed but rather that the input fit information is good or suspect. Exceptions would be more appropriate for failures to calculate the diagnostics.
The question introduced in this thread is designing a uniform interface for all diagnostic functions that is sufficiently robust for current (and immediate future diagnostics) that is compatible with the algorithm API.