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!

No-Recompile Mode

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!

VectorContinuousCallback

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 condition and affect! with structure -

using Pkg; Pkg.add("OrdinaryDiffEq")
using OrdinaryDiffEq
function condition(out,u,t,integrator)
  out[1] = u[1]
  out[2] = (u[3] - 10.0)u[3]
end

function affect!(integrator, idx)
  if idx == 1
    integrator.u[2] = -0.9integrator.u[2]
  elseif idx == 2
    integrator.u[4] = -0.9integrator.u[4]
  end
end

cb = VectorContinuousCallback(condition,affect!,2)

158.6s

Julia

VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(INITIALIZE_DEFAULT),Float64,Int64,Nothing}(condition, affect!, affect!, 2, INITIALIZE_DEFAULT, nothing, true, 10, Bool[true, true], 2.22045e-15, 0)

This is actually an example to simulate a ball bouncing between two walls and ground. We introduced out argument in condition and idx argument in affect!. What this does is, we make the user to set the output of the condition i in out[i]. Then whenever this becomes 0, we fire affect! with idx = i.

With the complete code we get -

function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -p
  du[3] = u[4]
  du[4] = 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")

2.7s

Julia

Next Plans

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!