In order to connect Stan with DifferentialEquations.jl the steps would be:
1. Create "diffeqcpp", a C++ interface to DifferentialEquations.jl (that would be similar to diffeqpy, diffeqr) possibly using CxxWrap.jl
2. Make it possible to evaluate vector-Jacobian products (VJP) with "diffeqcpp". Probably that would require ODE RHS to be coded as a string of Julia code, to make Julia AD libraries compatible with it.
At this point, it should be possible to call Julia solvers from C++ and evaluate the derivatives.
In Stan, there is stan::math::adj_jac_apply that makes it possible to define custom functions with custom VJP without having to deal with Stan autodiff types, it works for example with Eigen::Matrix<double>.
https://discourse.mc-stan.org/t/adj-jac-apply/5163
3. Make a class (let's call it JuliaODESolver) that implements two methods:
operator() // calls Julia solver for the given input
multiply_adjoint_jacobian() // evaluates VJP for the given vector
4. In .stan file add a custom function in "functions {}" block, and write a header file that implements that custom function. That would probably be one line
More info on using external C++ code is in Section 4.5 CmdStan Manual.
5. Modify cmdstan/main.cpp to initialize and finalize Julia context to be able to call Julia functions. This is probably the only place where Stan source itself needs to be modified.
I don't know what would be needed to make forward mode, and higher-order derivatives to work.
I think it would be much better for a fair benchmarking if there was a convenient and documented interface to Stan algorithms to use with user-provided log-density function, similar to DynamicHMC.jl and AdvancedHMC.jl libraries. It would be then easy to call it from Julia/Python/R/C++ or anything else.
1. Create "diffeqcpp", a C++ interface to DifferentialEquations.jl (that would be similar to diffeqpy, diffeqr) possibly using CxxWrap.jl
2. Make it possible to evaluate vector-Jacobian products (VJP) with "diffeqcpp". Probably that would require ODE RHS to be coded as a string of Julia code, to make Julia AD libraries compatible with it.
At this point, it should be possible to call Julia solvers from C++ and evaluate the derivatives.
In Stan, there is stan::math::adj_jac_apply that makes it possible to define custom functions with custom VJP without having to deal with Stan autodiff types, it works for example with Eigen::Matrix<double>. https://discourse.mc-stan.org/t/adj-jac-apply/5163
3. Make a class (let's call it JuliaODESolver) that implements two methods:
4. In .stan file add a custom function in "functions {}" block, and write a header file that implements that custom function. That would probably be one line
More info on using external C++ code is in Section 4.5 CmdStan Manual.
5. Modify cmdstan/main.cpp to initialize and finalize Julia context to be able to call Julia functions. This is probably the only place where Stan source itself needs to be modified.
I don't know what would be needed to make forward mode, and higher-order derivatives to work.
I think it would be much better for a fair benchmarking if there was a convenient and documented interface to Stan algorithms to use with user-provided log-density function, similar to DynamicHMC.jl and AdvancedHMC.jl libraries. It would be then easy to call it from Julia/Python/R/C++ or anything else.