Repository navigation
add functionality to interpolate in look-up tables #3365
Description
Activity
Interpolation has been proposed and discussed before: stan-dev/design-docs#18
I'm not sure why that work wasn't completed back then; I guess the author just moved on.Thanks! Well noted, and good to see I am not the first one to ask about it. Being new to Stan development, does this mean that the way forward is to try to reanimate or rewrite a design document for interpolation?
My main thoughts after browsing those ddiscussions are 1) rather than a new algorithm, I would wish to be able to use some subset of the most common methods available in other scientific computing libraries (matlab, scipy, boost, ...), and then 2) focus on fast evaluations.
I think the code I link to above is a step towards 1), but is not so impressive when it comes to speed. For one thing, it uses binary search to locate the interpolation interval, which one could beat if the input table has equally spaced data.
In that design doc, there is a comment
Ben Goodrich has pointed out that Boost is coming out with new interpolation schemes in an upcoming release.
It seems Boost has now these, and it would be great to use Boost (as we already use Boost for other things).
For naming I think it would be good to try to match the current style, and have for example
W = interpolate_1d_pchip_setup() interpolate_1d_pchip(W, ...)This name would be more clear than just
pchip, and would make it natural if other interpolation methods would be added.Alternatively, following Boost idea that many interpolators are splines, the names could be
W = interpolate_1d_pchip_setup() interpolate_1d_spline(W, ...)so that all setups that produce splines would use the same actual interpolation function. It would be fine to start implementing only pchip.
If using Boost, I think it would make sense to make a new design doc and link to the old one.
As you already have Stan implementation for PCHIP, you could also create a few tests that would pass or fail with certain inputs. This helps as then when making C++ functions, you need to add tests to Stan math side.
Reacted by Brian Ward and Martin LindénI agree with Aki on the naming and use of Boost. I actually think this probably falls below the threshold for what needs a design doc and would be happy to review a direct contribution, but more discussion is never unwelcome
Edit: It actually looks like the Boost implementation of PCHIP's precompute returns an opaque type rather than a matrix, so perhaps a reimplementation is better
Agree with all here, I think just a sketch of what the API will look like should be a way to start
@bmelinden do you need help with this? I can help
I can try to flesh out an API in some more detail. My understanding is also that Boost implements interpolators using a class structure, so that the pre-compute calls correspond to creators for those classes. In stan, I guess a pre-compute function that returns some existing Stan data type would be better. Hence, reimplementation.
For fast evaluation, cases with equally spaced grid points should be treated separately (faster evaluation). Boosts cardinal splines require equally spaced inputs. I think the output of the pre-calculation needs more than a matrix of polynomial coefficients: the grid points for non-equally spaced data, or step size and starting point for equally spaced data, so perhaps a tuple rather than a matrix?
The choice of interpolators in Boost makes a lot of sense to me, as does the naming pattern interpolate_1d_...
I think this argues for separate setup/eval pairs,
W = interpolate_1d_pchip_setup() yi=interpolate_1d_pchip(W, xi) V = interpolate_1d_cardinal_cubic_spline_setup() yj=interpolate_1d_cardinal_cubic_spline(V, xj)etc. Beyond this, I do need help finding my way around the codebase to get started adding tests etc, but I suspect OpenAI Codex could get me started on that.
Good points about using tuples and equally spaced grids. As we already use Boost integrators, AI agents should be able to help a lot in creating corresponding interpolators. I estimate I could add the interpolators with help from Claude in a few hours, but as we like to get new contributors I'm happy to let you try and if needed we can have a zoom call or you'll get feedback from @SteveBronder and @WardBrian for your PR
Sounds good! I started experimenting with OpenAI codex and pchip, and think I made some encouraging progress, so I would be happy to see if I can finish that and get back (might take a week or three).
Just letting you know I am still on it, after a false start (that maxed out my AI quota). I have a traft cubic hermite interpolator (since the boost pchip is a special case of hermite interpolation) and some tests (that all pass), but I have yet to test it out from stan. This is still mostly me learning the ropes of adding a new function and tests (and coding C++ with OpenRouter LLMs), the use case I think will be most useful are what boost calls ...cardinal... interpolators, which uses a constant spacing in the look-up table and therefore can be evaluated in constant time. Perhaps a good way forward would be to get the cubic hermite interpolator working as best I can, then collect some feedback on it, and then move on to add more interesting interpolators.
Good to hear! Happy to chat about this and look at any code you have so far. The main thing will be figuring out the reverse mode for propogating the adjoints from the interpolation. We will probably want something like we do for higher level functions like
integrate_1dor the ode solvers so we do not try to do autodiff through the interpolator itself
Many scientific and technical models incoporate tabulated data as part of the model definition. For example, to model temperature changes in some material where the heat capacity depends on temperature, one may need to compute the heat capacities at arbitrary temperatures even if measured values are only available at some specific temperatures. Many scientific computing tools implement various methods to interpolate in tabulated data, but Stan does not (as far as I know).
I think it would be good to add such functionality, as this would expand the scope of problems that can be modelled with Stan.
There are a lot of interpolation methods out there. To start with, I used OPenAI Codex to generate simple stan code for two common interpolation methods: PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) and cubic natural splines. Both have continuous first derivatives which I think is good for HMC simulations. (Piece-wise linear interpolation does not, and was not included). The code is available here: https://github.com/bmelinden/stan_interpolation, and include scripts that demonstrate both their use in Stan and good agreement with corresponding interpolation methods in R.
Both are implemented as two-step calls. First a setup using the tabulated data
W = pchip_setup(xk, yk);and then interpolated values can be calculated as
y = pchip_eval(x, W);Versions of these interpolation methods are available in for example Matlab (pchip, spline), Boost interpolation, scipy interpolate, and R (pracma pchip, splinefun). In fact, most software packages supply more interpolation methods than that, but these two may be a good enough start.
I would be happy to work more on this, but I am unsure on how to proceed.