48 ms·

Hi Chris, as a heads up, Stan actually has had three built-in methods for a while now. There is a non-stiff Adams-Moulton solver introduced in 2018. It unfortun

Hi Chris, as a heads up, Stan actually has had three built-in methods for a while now. There is a non-stiff Adams-Moulton solver introduced in 2018. It unfortunately was only just exposed in the Stan 2.23 documentation: https://mc-stan.org/docs/2_23/stan-users-guide/ode-solver-chapter.html https://mc-stan.org/docs/2_23/stan-users-guide/ode-solver-ch.... Certainly, the devs have been talking about adding more solvers for a while, including SDE and DDE solvers, and your DifferentialEquations.jl ecosystem is an excellent model; it is an area that we know Stan has been lacking in. I think Steve Gronder will be trying to work with you regarding benchmarking.

ChrisRackauckas 6y ago
Awesome. If we can get FFI with Stan I'd like to connect DifferentialEquations.jl to it and poke at it with various problems to see how well it does on a few things. We can provide custom gradients if there's an interface for it, but I couldn't figure out how to do it without modifying the Stan source itself.
- ivanyashchuk 6y ago
  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 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 return stan::math::adj_jac_apply<JuliaODESolver>(ode_solver_inputs); 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.

ChrisRackauckas 6y ago

Awesome. If we can get FFI with Stan I'd like to connect DifferentialEquations.jl to it and poke at it with various problems to see how well it does on a few things. We can provide custom gradients if there's an interface for it, but I couldn't figure out how to do it without modifying the Stan source itself.

ivanyashchuk 6y ago
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 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 return stan::math::adj_jac_apply<JuliaODESolver>(ode_solver_inputs); 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.