My GSoC with JuliaDiffEq - 1
Its really thrilling! I got a chance to spend my summer with the awesome community Julia is! Less Talk, More work: My Project is with JuliaDiffEq in improving the current benchmarks of the package in terms of performance. Basically it would be an impromptu project - with base requirements specified in the proposal.
During the community bonding period, I spent my most time reading Solving Ordinary Differential Equations - I (by Hairer, Ernst, Nørsett, Syvert P., Wanner, Gerhard) , developing the toolkit I would need for my project in the main project phase and finally fixing small issues including my pending pull requests. All of my contribution in this period can be divided into -
Developing a Benchmarking Server and Suite
We always have a problem in detecting regressions timely. Recently we were facing accuracy issues due to a particular commit, we all thought it was due to the refactoring of NLsolve routines but was actually due to a very small commit where we changed the initial value of
dt to handle
NaN in AD (automatic differentiation). To overcome such performance issues, it would be convenient to have a service which could tell us "Hey! This new cool feature is nice, but it causes these regressions". So this is what DiffEqBenchmarkServer is all about. I developed a neat server which is actually an abstraction over Jenkins, an open source CI/CD Pipeline. It works very similar to @JuliaRegistrator. All you have to do is to make a comment on the PR mentioning the bot @DiffEqBot to
runbenchmarks and DiffEqBot would do the work for you!
As all functions in Julia have unique type, parameterizing it causes many routines in the package to recompile everytime we change the function, even though the arguments and return types are same. So we have a neat package called FunctionalWrappers.jl which does the awesome work of wrapping the function in a structure correctly and then we pass it through our
solve, fixing this type instability. Unfortunately results were not that good for now (for instance, runtime regressions in stiff problems were prominent) so we have made it a keyword argument
norecompile and default its value to false. Once FunctionalWrappers becomes mature enough, we may change its default to true too!
DifferentialEquations already had a feature to plug user provided code when a particular event occurs, through callbacks. We also had
CallbackSet to chain up many callbacks together. But when you chain 1000+ callbacks together, it doesn't scale well enough. So we made a new type of callback called
VectorContinuousCallback where we have the function
affect! with structure -
using Pkg; Pkg.add("OrdinaryDiffEq") using OrdinaryDiffEq function condition(out,u,t,integrator) out = u out = (u - 10.0)u end function affect!(integrator, idx) if idx == 1 integrator.u = -0.9integrator.u elseif idx == 2 integrator.u = -0.9integrator.u end end cb = VectorContinuousCallback(condition,affect!,2)
This is actually an example to simulate a ball bouncing between two walls and ground. We introduced
out argument in
idx argument in
affect!. What this does is, we make the user to set the output of the condition
out[i]. Then whenever this becomes 0, we fire
idx = i.
With the complete code we get -
function f(du,u,p,t) du = u du = -p du = u du = 0.0 end u0 = [50.0,0.0,0.0,2.0] tspan = (0.0,15.0) p = 9.8 prob = ODEProblem(f,u0,tspan,p) sol = solve(prob,Tsit5(),callback=cb,dt=1e-3,adaptive=false) using Plots plt = plot(sol,vars=(3,1)) png(plt, "/results/ball.png")
By now, nothing from my proposal is done :p Everything I did in this month was discussed later and included many bug fixes. Now I aim to complete sub-projects of Implementing new time-step controllers and Pre-allocated Caches Feature before my next blog post. Till then, see ya!