C API refactor
I'd like to propose the following API for C (with help from @wardbrian and @roualdes). This will all be within the scope of extern "C"
.
Model and RNG holder
We encapsulate all the model-related information plus the random number generator state in a single class. It holds the RNG, but only a pointer to the model. The names and sizes will be instantiated lazily and can be freed along with the rest of the class.
struct model_rng {
stan::model::model_base* model_;
boost::ecuyer1988 prng_;
int param_unc_num_ = -1;
int param_num_ = -1;
int param_tp_num_ = -1;
int param_gq_num_ = -1;
int param_tp_gq_num_ = -1;
char* name_ = -1;
char* param_unc_names_ = nullptr;
char* param_names_ = nullptr;
char* param_tp_names_ = nullptr;
char* param_gq_names_ = nullptr;
char* param_tp_gq_names_ = nullptr;
};
The model will never be modified, but it cannot be declared const
because we don't have const correctness properly implemented on the Stan side. The PRNG can't be const
because it changes state. The other values can't be const
because they are lazily instantiated.
We need C functions to create and destroy the model from a data file. The construction function presupposes the model .hpp
file has been generated by stanc
and has been compiled using the CmdStan makefile---it will be linked in through the BridgeStan makefile.
model_rng* construct(char* data_file, unsigned int seed,
unsigned int chain_id);
void destruct(model_rng* mr);
Names
We have functions for the name of the model, and for the unconstrained and constrained parameters.
const char* name();
const char* param_names(model_rng* mr, bool include_tp, bool include_gq);
const char* param_unc_names(model_rng* mr);
The boolean flags indicate whether or not to include transformed parameters and generated quantities.
Parameters are printed as defined by the Stan model class. Parameter names are returned in CSV format without spaces. For example, the parameter theta[1, 2]
will be printed as theta.1.2
. Unconstrained parameters have either no indexes or have a single index as they are just used to map to the constrained parameters.
The const char*
that are returned are managed by the model_rng
class and its destruct
function.
Sizes
We have functions for the number of unconstrained parameters and number of constrained parameters with or without transformed parameters and generated quantities.
int param_num(model_rng* mr, bool include_tp, bool include_gq);
int param_unc_num(model_rng* mr);
Constraining and unconstraining transforms
For applied use, we will need unconstraining transforms to initialize from constrained parameter values, and we will need constraining transforms to produce the final results. The param_constrain
function has flags indicating whether to include the transformed parameters and/or generated quantities.
void param_constrain(model_rng* mr, bool include_tp, bool include_gq,
const double* theta_unc, double* theta);
void param_unconstrain(model_rng* mr, const double* theta, double* theta_unc);
void param_unconstrain_json(model_rng* mr, const char* json, double* theta_unc);
We might add yet another unconstrained method that supports Python dictionary specifications more naturally (though not with any Python dependencies in the C++). The signature is going to depend on the easiest way to get a Python dictionary to a Stan usable form that doesn't depend on reflecting an actual dictionary in C.
Log density, gradients, and Hessians
We need to be able to extract the density, gradient, and Hessian evaluated at an unconstrained parameter value. Once we've computed a gradient with Stan, the value is free. Similarly, once we've computed a Hessian, the gradient is free (other than the copy and alloc cost). The natural way to access seems to be with separate functions depending on which results we want.
double log_density(model_rng* mr, bool propto, bool jacobian,
const double* theta);
double log_density_gradient(model_rng* mr, bool propto, bool jacobian,
const double* theta, double* grad);
double log_density_hessian(model_rng* mr, bool propto, bool jacobian,
const double* theta, double* grad, double* hessian);
The functions return the log density directly rather than setting a pointer's value.
The Hessian method will use central finite diffs, which isn't very efficient or precise, but which works for any Stan model. We could alternatively swap in an autodiff implementation, but we'd need a compiler flag to signal that, because not every function supports autodiff Hessians. I probably won't get to this in the first refactor.