Solving Systems of Stochastic PDEs and using GPUs in Julia

What I want to describe in this post is how to solve stochastic PDEs in Julia using GPU parallelism. I will go from start to finish, describing how to use the type-genericness of the DifferentialEquations.jl library in order to write a code that uses within-method GPU-parallelism on the system of PDEs. This is mostly a proof of concept: the most efficient integrators for this problem are not compatible with GPU parallelism yet, and the GPU parallelism isn’t fully efficient yet. However, I thought it would be nice to show an early progress report showing that it works and what needs to be fixed in Base Julia and various libraries for us to get the full efficiency.

Our Problem: 2-dimensional Reaction-Diffusion Equations

The reaction-diffusion equation is a PDE commonly handled in systems biology which is a diffusion equation plus a nonlinear reaction term. The dynamics are defined as:

$u_t = D \Delta u + f(t,u)$

But this doesn’t need to only have a single “reactant” u: this can be a vector of reactants and the $f$ is then the nonlinear vector equations describing how these different pieces react together. Let’s settle on a specific equation to make this easier to explain. Let’s use a simple model of a 3-component system where A can diffuse through space to bind with the non-diffusive B to form the complex C (also non-diffusive, assume B is too big and gets stuck in a cell which causes C=A+B to be stuck as well). Other than the binding, we make each of these undergo a simple birth-death process, and we write down the equations which result from mass-action kinetics. If this all is meaningless to you, just understand that it gives the system of PDEs:

$A_t = D \Delta A + \alpha_A(x) - \beta_A A - r_1 A B + r_2 C$

$B_t = \alpha_B - \beta_B B - r_1 A B + r_2 C$

$C_t = \alpha_C - \beta_C C + r_1 A B - r_2 C$

One addition that was made to the model is that we let $\alpha_A(x)$ be the production of $A$, and we let that be a function of space so that way it only is produced on one side of our equation. Let’s make it a constant when x>80, and 0 otherwise, and let our spatial domain be $x \in [0,100]$ and $y \in [0,100]$.

This model is spatial: each reactant $u(t,x,y)$ is defined at each point in space, and all of the reactions are local, meaning that $f$ at spatial point $(x,y)$ only uses $u_i(t,x,y)$. This is an important fact which will come up later for parallelization.

Discretizing the PDE into ODEs

In order to solve this via a method of lines (MOL) approach, we need to discretize the PDE into a system of ODEs. Let’s do a simple uniformly-spaced grid finite difference discretization. Choose $dx = 1$ and $dy = 1$ so that we have 100*100=10000 points for each reactant. Notice how fast that grows! Put the reactants in a matrix such that A[i,j] = $A(x_j,y_i)$, i.e. the columns of the matrix is the $x$ values and the rows are the $y$ values (this way looking at the matrix is essentially like looking at the discretized space).

So now we have 3 matrices (A, B, and C) for our reactants. How do we discretize the PDE? In this case, the diffusion term simply becomes a tridiagonal matrix $M$ where $[1,-2,1]$ is central band. You can notice that $MA$ performs diffusion along the columns of $A$, and so this is diffusion along the $y$. Similarly, $AM$ flips the indices and thus does diffusion along the rows of $A$ making this diffusion along $x$. Thus $D(M_yA + AM_x)$ is the discretized Laplacian (we could have separate diffusion constants and $dx \neq dy$ if we want by using different constants on the $M$, but let’s not do that for this simple example. I’ll leave that as an exercise for the reader). I enforced a Neumann boundary condition with zero derivative (also known as a no-flux boundary condition) by reflecting the changes over the boundary. Thus the derivative operator is generated as:

const Mx = full(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
const My = copy(Mx)
# Do the reflections, different for x and y operators
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

I also could have done this using the DiffEqOperators.jl library, but I wanted to show what it truly is at its core.

Since all of the reactions are local, we only have each point in space react separately. Thus this represents itself as element-wise equations on the reactants. Thus we can write it out quite simply. The ODE which then represents the PDE is thus in pseudo Julia code:

DA = D*(M*A + A*M)
@. DA + α₁ - β₁*A - r₁*A*B + r₂*C
@. α₂ - β₂*B - r₁*A*B + r₂*C
@. α₃ - β₃*C + r₁*A*B - r₂*C

Note here that I am using α₁ as a matrix (or row-vector, since that will broadcast just fine) where every point in space with x<80 has this zero, and all of the others have it as a constant. The other coefficients are all scalars. How do we do this with the ODE solver?

Our Type: ArrayPartition

The ArrayPartition is an interesting type from RecursiveArrayTools.jl which allows you to define “an array” as actually being different discrete subunits of arrays. Let’s assume that our initial condition is zero for everything and let the production terms build it up. This means that we can define:

A = zeros(M,N); B  = zeros(M,N); C = zeros(M,N)

Now we can put them together as:

u0 = ArrayPartition((A,B,C))

You can read the RecursiveArrayTools.jl README to get more familiar with what the ArrayPartition is, but really it’s an array where u[i] indexes into A first, B second, then C. It also has efficient broadcast, doing the A, B and C parts together (and this is efficient even if they don’t match types!). But since this acts as an array, to DifferentialEquations.jl it is an array!

The important part is that we can “decouple” the pieces of the array at anytime by accessing u.x, which holds our tuple of arrays. Thus our ODE using this ArrayPartition as its container can be written as follows:

function f(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
DA = D*(M*A + A*M)
@. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

where this is using @. to do inplace updates on our du to say how the full ArrayPartition should update in time. Note that we can make this more efficient by adding some cache variables to the diffusion matrix multiplications and using A_mul_B!, but let’s ignore that for now.

Together, the ODE which defines our PDE is thus:

prob = ODEProblem(f,u0,(0.0,100.0))
sol = solve(prob,BS3())

if I want to solve it on $t \in [0,100]$. Done! The solution gives back ArrayPartitions (and interpolates to create new ones if you use sol(t)). We can plot it in Plots.jl

and see the pretty gradients. Using this 3rd order explicit adaptive Runge-Kutta method we solve this equation in about 40 seconds. That’s okay.

Some Optimizations

There are some optimizations that can still be done. When we do A*B as matrix multiplication, we create another temporary matrix. These allocations can bog down the system. Instead we can pre-allocate the outputs and use the inplace functions A_mul_B! to make better use of memory. The easiest way to store these cache arrays are constant globals, but you can use closures (anonymous functions which capture data, i.e. (x)->f(x,y)) or call-overloaded types to do it without globals. The globals way (the easy way) is simply:

const MyA = zeros(N,N)
const AMx = zeros(N,N)
const DA = zeros(N,N)
function f(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
A_mul_B!(MyA,My,A)
A_mul_B!(AMx,A,Mx)
@. DA = D*(MyA + AMx)
@. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

For reference, closures looks like:

MyA = zeros(N,N)
AMx = zeros(N,N)
DA = zeros(N,N)
function f_full(t,u,du,MyA,AMx,DA)
A,B,C = u.x
dA,dB,dC = du.x
A_mul_B!(MyA,My,A)
A_mul_B!(AMx,A,Mx)
@. DA = D*(MyA + AMx)
@. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end
f = (t,u,du)-> f_full(t,u,du,MyA,AMx,DA)

and a call overloaded type looks like:

struct MyFunction{T} <: Function
MyA::T
AMx::T
DA::T
end

function (ff::MyFunction)(t,u,du)
# This is a function which references itself via ff
A,B,C = u.x
dA,dB,dC = du.x
A_mul_B!(ff.MyA,My,A)
A_mul_B!(ff.AMx,A,Mx)
@. ff.DA = D*(ff.MyA + ff.AMx)
@. dA = f.DA + α₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

MyA = zeros(N,N)
AMx = zeros(N,N)
DA = zeros(N,N)

f = MyFunction(MyA,AMx,DA)
# Now f(t,u,du) is our function!

These last two ways enclose the pointer to our cache arrays locally but still present a function f(t,u,du) to the ODE solver.

Now since PDEs are large, many times we don’t care about getting the whole timeseries. Using the output controls from DifferentialEquations.jl, we can make it only output the final timepoint.

sol = solve(prob,BS3(),progress=true,save_everystep=false,save_start=false)

Also, if you’re using Juno this’ll give you a nice progress bar so you can track how it’s going.

We are using an explicit Runge-Kutta method here because that’s what works with GPUs so far. Matrix factorizations need to be implemented for GPUArrays before the implicit (stiff) solvers will be available, so here we choose BS3 since it’s fully broadcasting (not all methods are yet) and it’s fully GPU compatible. In practice, right now using an NxNx3 tensor as the initial condition / dependent variable with either OrdinaryDiffEq’s Rosenbrock23(), Rodas4(), or Sundials’ CVODE_BDF() is actually more efficient right now. But after Julia fixes its broadcasting issue and with some updates to Julia’s differentiation libraries to handle abstract arrays like in DiffEqDiffTools.jl, the stiff solvers will be usable with GPUs and all will be well.

Thus for reference I will show some ways to do this efficiently with stiff solvers. With a stiff solver we will not want to factorize the dense Jacobian since that would take forever. Instead we can use something like Sundials’ Krylov method:

u0 = zeros(N,N,3)
const MyA = zeros(N,N);
const AMx = zeros(N,N);
const DA = zeros(N,N)
function f(t,u,du)
A = @view u[:,:,1]
B = @view u[:,:,2]
C = @view u[:,:,3]
dA = @view du[:,:,1]
dB = @view du[:,:,2]
dC = @view du[:,:,3]
A_mul_B!(MyA,My,A)
A_mul_B!(AMx,A,Mx)
@. DA = D*(MyA + AMx)
@. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

# Solve the ODE
prob = ODEProblem(f,u0,(0.0,100.0))
using Sundials
@time sol = solve(prob,CVODE_BDF(linear_solver=:BCG))

and that will solve it in about a second. In this case it wouldn’t be more efficient to use the banded linear solver since the system of equations tends to have different parts of the system interact which makes the bands large, and thus a Krylov method is preferred. See this part of the docs for details on the available linear solvers from Sundials. DifferentialEquations.jl exposes a ton of Sundials’ possible choices so hopefully one works for your problem (preconditioners coming soon).

To do something similar with OrdinaryDiffEq.jl, we would need to make use of the linear solver choices in order to override the internal linear solve functions with some kind of sparse matrix solver like a Krylov method from IterativeSolvers.jl. For this size of problem though a multistep method like BDF is probably preferred though, at least until we implement some IMEX methods.

So if you want to solve it quickly right now, that’s how you do it. But let’s get back to our other story: the future is more exciting.

The Full ODE Code

As a summary, here’s a full PDE code:

using OrdinaryDiffEq, RecursiveArrayTools

# Define the constants for the PDE
const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const D = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 100
const X = reshape([i for i in 1:100 for j in 1:100],N,N)
const Y = reshape([j for i in 1:100 for j in 1:100],N,N)
const α₁ = 1.0.*(X.>=80)

const Mx = full(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

# Define the initial condition as normal arrays
A = zeros(N,N); B  = zeros(N,N); C = zeros(N,N)
u0 = ArrayPartition((A,B,C))

const MyA = zeros(N,N);
const AMx = zeros(N,N);
const DA = zeros(N,N)
# Define the discretized PDE as an ODE function
function f(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
A_mul_B!(MyA,My,A)
A_mul_B!(AMx,A,Mx)
@. DA = D*(MyA + AMx)
@. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

# Solve the ODE
prob = ODEProblem(f,u0,(0.0,100.0))
sol = solve(prob,BS3(),progress=true,save_everystep=false,save_start=false)

using Plots; pyplot()
p1 = surface(X,Y,sol[end].x[1],title = "[A]")
p2 = surface(X,Y,sol[end].x[2],title = "[B]")
p3 = surface(X,Y,sol[end].x[3],title = "[C]")
plot(p1,p2,p3,layout=grid(3,1))

Making Use of GPU Parallelism

That was all using the CPU. How do we make turn on GPU parallelism with DifferentialEquations.jl? Well, you don’t. DifferentialEquations.jl “doesn’t have GPU bits”. So wait… can we not do GPU parallelism? No, this is the glory of type-genericness, especially in broadcasted operations. To make things use the GPU, we simply use a GPUArray. If instead of zeros(N,M) we used GPUArray(zeros(N,M)), then u becomes an ArrayPartition of GPUArrays. GPUArrays naturally override broadcast such that dotted operations are performed on the GPU. DifferentialEquations.jl uses broadcast internally (except in this list of current exceptions due to a limitation with Julia’s inference engine which I have discussed with Jameson Nash (@vtjnash) who mentioned this should be fixed in Julia’s 1.0 release), and thus just by putting the array as a GPUArray, the array-type will take over how all internal updates are performed and turn this algorithm into a fully GPU-parallelized algorithm that doesn’t require copying to the CPU. Wasn’t that simple?

From that you can probably also see how to multithread everything, or how to set everything up with distributed parallelism. You can make the ODE solvers do whatever you want by defining an array type where the broadcast does whatever special behavior you want.

So to recap, the entire difference from above is changing to:

using CLArrays
gA = CLArray(A); gB  = CLArray(B); gC = CLArray(C)
const gMx = CLArray(Mx)
const gMy = CLArray(My)
const gα₁ = CLArray(α₁)
gu0 = ArrayPartition((gA,gB,gC))

const gMyA = zeros(N,N)
const gAMx = zeros(N,N)
const gDA = zeros(N,N)
function gf(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
A_mul_B!(gMyA,gMy,A)
A_mul_B!(gAMx,A,gMx)
@. DA = D*(gMyA + AgMx)
@. dA = DA + gα₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

prob2 = ODEProblem(gf,gu0,(0.0,100.0))
GPUArrays.allowslow(false) # makes sure none of the slow fallbacks are used

prob2 = ODEProblem(gf,gu0,(0.0,100.0))
sol = solve(prob2,BS3(),progress=true,save_everystep=false,save_start=false)
# Adaptivity currently fails due to https://github.com/JuliaGPU/CLArrays.jl/issues/10

You can use CUArrays if you want as well. It looks exactly the same as using CLArrays except you exchange the CLArray calls to CUArray. Go have fun.

And Stochastic PDEs?

Why not make it an SPDE? All that we need to do is extend each of the PDE equations to have a noise function. In this case, let’s use multiplicative noise on each reactant. This means that our noise update equation is:

function g(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
@. dA = γ₁*A
@. dB = γ₂*A
@. dC = γ₃*A
end

Now we just define and solve the system of SDEs:

prob = SDEProblem(f,g,u0,(0.0,100.0))
sol = solve(prob,SRIW1())

We can see the cool effect that diffusion dampens the noise in [A] but is unable to dampen the noise in [B] which results in a very noisy [C]. The stiff SPDE takes much longer to solve even using high order plus adaptivity because stochastic problems are just that much more difficult (current research topic is to make new algorithms for this!). It gets GPU’d just by using GPUArrays like before. But there we go: solving systems of stochastic PDEs using high order adaptive algorithms with within-method GPU parallelism. That’s gotta be a first? The cool thing is that nobody ever had to implement the GPU-parallelism either, it just exists by virtue of the Julia type system.

Side Notes

Warning: This can take awhile to solve! An explicit Runge-Kutta algorithm isn’t necessarily great here, though to use a stiff solver on a problem of this size requires once again smartly choosing sparse linear solvers. The high order adaptive method is pretty much necessary though since something like Euler-Maruyama is simply not stable enough to solve this at a reasonable dt. Also, the current algorithms are not so great at handling this problem. Good thing there’s a publication coming along with some new stuff…

Note: the version of SRIW1 which uses broadcast for GPUs is not on the current versions of StochasticDiffEq.jl since it’s slower due to a bug when fusing too many broadcasts which will hopefully get fixed in one of Julia’s 1.x releases. Until then, GPUs cannot be used with this algorithm without a (quick) modification.

Conclusion

So that’s where we’re at. GPU parallelism works because of abstract typing. But in some cases we need to help the GPU array libraries get up to snuff to handle all of the operations, and then we’ll really be in business! Of course there’s more optimizing that needs to be done, and we can do this by specializing code paths on bottlenecks as needed.

I think this is at least a nice proof of concept showing that Julia’s generic algorithms allow for one to not only take advantage of things like higher precision, but also take advantage of parallelism and extra hardware without having to re-write the underlying algorithm. There’s definitely more work that needs to be done, but I can see this usage of abstract array typing as being one of Julia’s “killer features” in the coming years as the GPU community refines its tools. I’d give at least a year before all of this GPU stuff is compatible with stiff solvers and linear solver choices (so that way it can make use of GPU-based Jacobian factorizations and Krylov methods). And comparable methods for SDEs are something I hope to publish soon since the current tools are simply not fit for this scale of problem: high order, adaptivity, sparse linear solvers, and A/L-stability all need to be combined in order to tackle this problem efficiently.

Full Script

Here’s the full script for recreating everything:

#######################################################
### Solve the PDE
#######################################################

using OrdinaryDiffEq, RecursiveArrayTools

# Define the constants for the PDE
const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const D = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 100
const X = reshape([i for i in 1:100 for j in 1:100],N,N)
const Y = reshape([j for i in 1:100 for j in 1:100],N,N)
const α₁ = 1.0.*(X.>=80)

const Mx = full(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

# Define the initial condition as normal arrays
A = zeros(N,N); B  = zeros(N,N); C = zeros(N,N)
u0 = ArrayPartition((A,B,C))

const MyA = zeros(N,N);
const AMx = zeros(N,N);
const DA = zeros(N,N)
# Define the discretized PDE as an ODE function
function f(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
A_mul_B!(MyA,My,A)
A_mul_B!(AMx,A,Mx)
@. DA = D*(MyA + AMx)
@. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

# Solve the ODE
prob = ODEProblem(f,u0,(0.0,100.0))
@time sol = solve(prob,BS3(),progress=true,save_everystep=false,save_start=false)

using Plots; pyplot()
p1 = surface(X,Y,sol[end].x[1],title = "[A]")
p2 = surface(X,Y,sol[end].x[2],title = "[B]")
p3 = surface(X,Y,sol[end].x[3],title = "[C]")
plot(p1,p2,p3,layout=grid(3,1))

#######################################################
### Solve the PDE using CLArrays
#######################################################

using CLArrays
gA = CLArray(A); gB  = CLArray(B); gC = CLArray(C)
const gMx = CLArray(Mx)
const gMy = CLArray(My)
const gα₁ = CLArray(α₁)
gu0 = ArrayPartition((gA,gB,gC))

const gMyA = CLArray(MyA)
const gAMx = CLArray(AMx)
const gDA = CLArray(DA)
function gf(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
A_mul_B!(gMyA,gMy,A)
A_mul_B!(gAMx,A,gMx)
@. gDA = D*(gMyA + gAMx)
@. dA = gDA + gα₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

prob2 = ODEProblem(gf,gu0,(0.0,100.0))
GPUArrays.allowslow(false)

prob2 = ODEProblem(gf,gu0,(0.0,100.0))
sol = solve(prob2,BS3(),progress=true,save_everystep=false,save_start=false)
# Adaptivity currently fails due to https://github.com/JuliaGPU/CLArrays.jl/issues/10

#######################################################
### Solve the SPDE
#######################################################

using StochasticDiffEq

function g(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
@. dA = γ₁*A
@. dB = γ₂*A
@. dC = γ₃*A
end

prob3 = SDEProblem(f,g,u0,(0.0,100.0))
sol = solve(prob3,SRIW1(),progress=true,save_everystep=false,save_start=false)

p1 = surface(X,Y,sol[end].x[1],title = "[A]")
p2 = surface(X,Y,sol[end].x[2],title = "[B]")
p3 = surface(X,Y,sol[end].x[3],title = "[C]")
plot(p1,p2,p3,layout=grid(3,1))

# Exercise: Do SPDE + GPU

The post Solving Systems of Stochastic PDEs and using GPUs in Julia appeared first on Stochastic Lifestyle.

Data Wrangling in Julia based on dplyr Flights Tutorials

A couple of my favorite tutorials for wrangling data in R with dplyr are Hadley Wickham’s dplyr package vignette and Kevin Markham’s dplyr tutorial. I enjoy the tutorials because they concisely illustrate how to use a small set of verb-based functions to carry out common data wrangling tasks.

I tend to use Python to wrangle data, but I’m exploring the Julia programming language so I thought creating a similar dplyr-based tutorial in Julia would be a fun way to examine Julia’s capabilities and syntax. Julia has several packages that make it easier to deal with tabular data, including DataFrames and DataFramesMeta.

The DataFrames package provides functions for reading and writing, split-apply-combining, reshaping, joining, sorting, querying, and grouping tabular data. The DataFramesMeta package provides a set of macros that are similar to dplyr’s verb-based functions in that they offer a more convenient, readable syntax for munging data and chaining together multiple operations.

Data

For this tutorial, let’s following along with Kevin’s tutorial and use the hflights dataset. You can obtain the dataset from R with the following commands or simply download it here: hflights.csv

install.packages("hflights") library(hflights) write.csv(hflights, "hflights.csv")

To begin, let’s start the Julia REPL, load the DataFrames and DataFramesMeta packages, and load and inspect the hflights dataset:

using DataFrames using DataFramesMeta

hflights = readtable("/Users/clinton/Documents/Julia/hflights.csv"); size(hflights) names(hflights) head(hflights) describe(hflights)

The semicolon on the end of the readtable command prevents it from printing the dataset to the screen. The size command returns the number of rows and columns in the dataset. You can specify you only want the number of rows with size(hflights, 1) or columns with size(hflights, 2). This dataset contains 227,496 rows and 21 columns. The names command lists the column headings. By default, the head command prints the header row and six data rows. You can specify the number of data rows to display by adding a second argument, e.g. head(hflights, 10). The describe command prints summary statistics for each column.

@where: Keep rows matching criteria

AND: All of the conditions must be true for the returned rows

# Julia DataFrames approach to view all flights on January 1 hflights[.&(hflights[:Month] .== 1, hflights[:DayofMonth] .== 1), :]

# DataFramesMeta approach @where(hflights, :Month .== 1, :DayofMonth .== 1)

Julia’s DataFrames’ row filtering syntax is similar to R’s syntax. To specify multiple AND conditions, use “.&()” and place the filtering conditions, separated by commas, between the parentheses. Like dplyr’s filter function, DataFramesMeta’s @where macro simplifies the syntax and makes the command easier to read.

OR: One of the conditions must be true for the returned rows

# Julia DataFrames approach to view all flights where either AA or UA is the carrier hflights[.|(hflights[:UniqueCarrier] .== "AA", hflights[:UniqueCarrier] .== "UA"), :]

# DataFramesMeta approach @where(hflights, .|(:UniqueCarrier .== "AA", :UniqueCarrier .== "UA"))

To specify multiple OR conditions, use “.|()” and place the filtering conditions between the parentheses. Again, the DataFramesMeta approach is more concise.

SET: The values in a column are in a set of interest

# Julia DataFrames approach to view all flights where the carrier is in Set(["AA", "UA"]) carriers_set = Set(["AA", "UA"]) hflights[findin(hflights[:UniqueCarrier], carriers_set), :]

# DataFramesMeta approach @where(hflights, findin(:UniqueCarrier, carriers_set))

To filter for rows where the values in a particular column are in a specific set of interest, create a Set with the values you’re interested in and then specify the column and your set of interest in the findin function.

PATTERN / REGULAR EXPRESSION: The values in a column match a pattern

# Julia DataFrames approach to view all flights where the carrier matches the regular expression r"AA|UA" carriers_pattern = r"AA|UA" hflights[[ismatch(carriers_pattern, String(carrier)) for carrier in hflights[:UniqueCarrier]], :]

# DataFramesMeta approach @where(hflights, [ismatch(carriers_pattern, String(carrier)) for carrier in :UniqueCarrier])

To filter for rows where the values in a particular column match a pattern, create a regular expression and then use it in the ismatch function in an array comprehension.

@select: Pick columns by name

# Julia DataFrames approach to selecting columns hflights[:, [:DepTime, :ArrTime, :FlightNum]]

# DataFramesMeta approach @select(hflights, :DepTime, :ArrTime, :FlightNum)

Julia’s DataFrames’ syntax for selecting columns is similar to R’s syntax. Like dplyr’s select function, DataFramesMeta’s @select macro simplifies the syntax and makes the command easier to read.

# Julia DataFrames approach to selecting columns # first three columns hflights[:, 1:3] # pattern / regular expression heading_pattern = r"Taxi|Delay" hflights[:, [ismatch(heading_pattern, String(name)) for name in names(hflights)]] # startswith hflights[:, filter(name -> startswith(String(name), "Arr"), names(hflights))] # endswith hflights[:, filter(name -> endswith(String(name), "Delay"), names(hflights))] # contains hflights[:, filter(name -> contains(String(name), "Month"), names(hflights))]

# AND conditions hflights[:, filter(name -> startswith(String(name), "Arr") && endswith(String(name), "Delay"), names(hflights))] # OR conditions hflights[:, filter(name -> startswith(String(name), "Arr") || contains(String(name), "Cancel"), names(hflights))]

# DataFramesMeta approach # first three columns @select(hflights, 1:3) # pattern / regular expression heading_pattern = r"Taxi|Delay" @select(hflights, [ismatch(heading_pattern, String(name)) for name in names(hflights)]) # startswith @select(hflights, filter(name -> startswith(String(name), "Arr"), names(hflights))) # endswith @select(hflights, filter(name -> endswith(String(name), "Delay"), names(hflights))) # contains @select(hflights, filter(name -> contains(String(name), "Month"), names(hflights)))

# AND conditions @select(hflights, filter(name -> startswith(String(name), "Arr") && endswith(String(name), "Delay"), names(hflights))) # OR conditions @select(hflights, filter(name -> startswith(String(name), "Arr") || contains(String(name), "Cancel"), names(hflights)))

# Kevin Markham's multiple select conditions example # select(flights, Year:DayofMonth, contains("Taxi"), contains("Delay")) # Julia Version of Kevin's Example # Taxi or Delay in column heading mask = [ismatch(r"Taxi|Delay", String(name)) for name in names(hflights)] # Also include first three columns, i.e. Year, Month, DayofMonth mask[1:3] = true @select(hflights, mask)

These examples show you can select columns by position and name, and you can combine multiple conditions with AND, “&&”, or OR, “||”. Similar to filtering rows, you can select specific columns based on a pattern by using the ismatch function in an array comprehension. You can also use contains, startswith, and endswith in the filter function to select columns that contain, start with, or end with a specific text pattern.

“Chaining” or “Pipelining”

In R, dplyr provides, via the magrittr package, the %>% operator, which enables you to chain together multiple commands into a single data transformation pipeline in a very readable fashion. In Julia, the DataFramesMeta package provides the @linq macro and |> symbol to enable similar functionality. Alternatively, you can load the Lazy package and use an @> begin end block to chain together multiple commands.

# Chaining commands with DataFrameMeta’s @linq macro @linq hflights[find(.!isna.(hflights[:,:DepDelay])), :] |> @where(:DepDelay .> 60) |> @select(:UniqueCarrier, :DepDelay)

# Chaining commands with Lazy’s @> begin end block using Lazy @> begin hflights[find(.!isna.(hflights[:,:DepDelay])), :] @where(:DepDelay .> 60) @select(:UniqueCarrier, :DepDelay) end

These two blocks of code produce the same result, a DataFrame containing carrier names and departure delays for which the departure delay is greater than 60. In each chain, the first expression is the input DataFrame, e.g. hflights. In these examples, I use the find and !isna. functions to start with a DataFrame that doesn’t contain NA values in the DepDelay column because the commands fail when NAs are present. I prefer the @linq macro version over the @> begin end version because it’s so similar to the dplyr-magrittr syntax, but both versions are more succinct and readable than their non-chained versions. The screen shot shows how to assign the pipeline results to variables.

@orderby: Reorder rows

Both DataFrames and DataFramesMeta provide functions for sorting rows in a DataFrame by values in one or more columns. In the first pair of examples, we want to select the UniqueCarrier and DepDelay columns and then sort the results by the values in the DepDelay column in descending order. The last example shows how to sort by multiple columns with the @orderby macro.

# Julia DataFrames approach to sorting sort(hflights[find(.!isna.(hflights[:,:DepDelay])), [:UniqueCarrier, :DepDelay]], cols=[order(:DepDelay, rev=true)])

# DataFramesMeta approach (add a minus sign before the column symbol for descending) @linq hflights[find(.!isna.(hflights[:,:DepDelay])), :] |> @select(:UniqueCarrier, :DepDelay) |> @orderby(-:DepDelay)

# Sort hflights dataset by Month, descending, and then by DepDelay, ascending @linq hflights |> @orderby(-:Month, :DepDelay)

DataFrames provides the sort and sort! functions for ordering rows in a DataFrame. sort! orders the rows, inplace. The DataFrames user guide provides additional examples of ordering rows, in ascending and descending order, based on multiple columns, as well as applying functions to columns, e.g. uppercase, before using the column for sorting.

DataFramesMeta provides the @orderby macro for ordering rows in a DataFrame. Specify multiple column names in the @orderby macro to sort the rows by multiple columns. Use a minus sign before a column name to sort in descending order.

Creating new variables in Julia DataFrames is similar to creating new variables in Python and R. You specify a new column name in square brackets after the name of the DataFrame and assign it a collection of values, sometimes based on values in other columns. DataFramesMeta’s @transform macro simplifies the syntax and makes the transformation more readable.

# Julia DataFrames approach to creating new variable hflights[:Speed] = hflights[:Distance] ./ hflights[:AirTime] .* 60 hflights[:, [:Distance, :AirTime, :Speed]]

# Delete the variable so we can recreate it with DataFramesMeta approach delete!(hflights, :Speed)

# DataFramesMeta approach @linq hflights |> @select(:Distance, :AirTime) |> @transform(Speed = :Distance ./ :AirTime .* 60) |> @select(:Distance, :AirTime, :Speed)

# Save the new column in the original DataFrame hflights = @linq hflights |> @transform(Speed = :Distance ./ :AirTime .* 60)

The first code block illustrates how to create a new column in a DataFrame and assign it values based on values in other columns. The second code block shows you can use delete! to delete a column. The third example demonstrates the DataFramesMeta approach to creating a new column using the @transform macro. The last example shows how to save a new column in an existing DataFrame using the @transform macro by assigning the result of the transformation to the existing DataFrame.

@by: Reduce variables to values (Grouping and Summarizing)

dplyr provides group_by and summarise functions for grouping and summarising data. DataFrames and DataFramesMeta also support the split-apply-combine strategy with the by function and the @by macro, respectively. Here Julia versions of Kevin’s summarise examples.

# Julia DataFrames approach to grouping and summarizing by(hflights[complete_cases(hflights[[Symbol(name) for name in names(hflights)]]), :], :Dest, df -> DataFrame(meanArrDelay = mean(df[:ArrDelay])))

# DataFramesMeta approach @linq hflights[complete_cases(hflights[[Symbol(name) for name in names(hflights)]]), :] |> @by(:Dest, meanArrDelay = mean(:ArrDelay))

DataFrames and DataFramesMeta don’t have dplyr’s summarise_each function, but it’s easy to apply different functions to multiple columns inside the @by macro.

@linq hflights |> @by(:UniqueCarrier, meanCancelled = mean(:Cancelled), meanDiverted = mean(:Diverted))

@linq hflights[complete_cases(hflights[[Symbol(name) for name in names(hflights)]]), :] |> @by(:UniqueCarrier, minArrDelay = minimum(:ArrDelay), maxArrDelay = maximum(:ArrDelay), minDepDelay = minimum(:DepDelay), maxDepDelay = maximum(:DepDelay))

DataFrames and DataFramesMeta also don’t have dplyr’s n and n_distinct functions, but you can count the number of rows in a group with size(df, 1) or nrow(df), and you can count the number of distinct values in a group with countmap.

# Group by Month and DayofMonth, count the number of flights, and sort descending # Count the number of rows with size(df, 1) sort(by(hflights, [:Month,:DayofMonth], df -> DataFrame(flight_count = size(df, 1))), cols=[order(:flight_count, rev=true)])

# Group by Month and DayofMonth, count the number of flights, and sort descending # Count the number of rows with nrow(df) sort(by(hflights, [:Month,:DayofMonth], df -> DataFrame(flight_count = nrow(df))), cols=[order(:flight_count, rev=true)])

# Split grouping and sorting into two separate operations g = by(hflights, [:Month,:DayofMonth], df -> DataFrame(flight_count = nrow(df))) sort(g, cols=[order(:flight_count, rev=true)])

# For each destination, count the total number of flights and the number of distinct planes by(hflights[find(.!isna.(hflights[:,:TailNum])),:], :Dest) do df DataFrame(flight_count = size(df,1), plane_count = length(keys(countmap(df[:,:TailNum])))) end

While these examples reproduce the results in Kevin’s dplyr tutorial, they’re definitely not as succinct and readable as the dplyr versions. Grouping by multiple columns, summarizing with counts and distinct counts, and gracefully chaining these operations are areas where DataFrames and DataFramesMeta can improve.

Other useful convenience functions

Randomly sampling a fixed number or fraction of rows from a DataFrame can be a helpful operation. dplyr offers the sample_n and sample_frac functions to perform these operations. In Julia, StatsBase provides the sample function, which you can repurpose to achieve similar results.

 using StatsBase # randomly sample a fixed number of rows hflights[sample(1:nrow(hflights), 5), :] hflights[sample(1:size(hflights,1), 5), :]

# randomly sample a fraction of rows hflights[sample(1:nrow(hflights), ceil(Int,0.0001*nrow(hflights))), :] hflights[sample(1:size(hflights,1), ceil(Int,0.0001*size(hflights,1))), :]

Randomly sampling a fixed number of rows is fairly straightforward. You use the sample function to randomly select a fixed number of rows, in this case five, from the DataFrame. Randomly sampling a fraction of rows is slightly more complicated because, since the sample function takes an integer for the number of rows to return, you need to use the ceil function to convert the fraction of rows, in this case 0.0001*nrow(hflights), into an integer.

Conclusion

In R, dplyr sets a high bar for wrangling data well with succinct, readable code. In Julia, DataFrames and DataFramesMeta provide many useful functions and macros that produce similar results; however, some of the syntax isn’t as concise and clear as it is with dplyr, e.g. selecting columns in different ways and chaining together grouping and summarizing operations. These are areas where Julia’s packages can improve.

I enjoyed becoming more familiar with Julia by reproducing much of Kevin’s dplyr tutorial. It was also informative to see differences in functionality and readability between dplyr and Julia’s packages. I hope you enjoyed this tutorial and find it to be a useful reference for wrangling data in Julia.

Filed under: Analytics, General, Julia, Python, R, Statistics Tagged: DataFrames, DataFramesMeta, dplyr, Julia, Python, R

DifferentialEquations.jl 3.0 and a Roadmap for 4.0

I am pleased to announce the release of DifferentialEquations.jl 3.0. In the last DiffEq blog post I described the current state of JuliaDiffEq and DifferentialEquations.jl along with the trajectory that we hoped to take. We identified (at that time) current shortcomings of the software and our plans to remedy them. I also recently did a survey of differential equation suites in order to understand where we stand and see where we need to improve. These research efforts were used to put together a list of goals that were systematically achieved during 3.0. What I would like to do this time around is give a broad overview of what we have released in the 3.0 timeframe, the goals that we have achieved, and the goals that we are putting off (for next Google Summer of Code?). And then, more importantly, I want to set some milestones for the next version. If you want to dig into our new features and start using them, please see the documentation. If you want to read the release posts, see the official JuliaDiffEq blog

A Quick Review of DifferentialEquations.jl Pre-3.0

In 1.0, we made every thing work with generic types and event handling. In all of the native Julia solvers you could use arbitrary arithmetic and use events to have the ODEs do crazy things like change size over time. This was about features. In 2.0, we expanded our capabilities to cover “most” of what users tend to need. A broad array of ordinary differential equation (ODE) solvers, a broad array of stochastic differential equation (SDE) solvers, delay differential equation (DDE) solvers, and some partial differential equation (PDE) solvers. We added addons for parameter estimation, sensitivity analysis, uncertainty quantification, etc. This was really exciting because it was the first set of differential equation solvers which had this range of applicability. What this did was make it possible to solve many different types of problems. You “could” solve the problems. There were some edge cases for sure, but the main areas where the vast majority of individuals were looking was hit. But there were two major remaining warts: stiff problems and PDEs.

Introducing DifferentialEquations.jl 3.0

There was an issue. There are some specific types of problems, namely stiff differential equations, which require specific types of methods. We had wrappers to common C/Fortran solver for these, but this meant that we lost the type flexibility and event handling when solving these equations. We couldn’t handle some of the more difficult problems like state-dependent delays as well. Thus these types of problems were the focus of 3.0: to have some semblance of “completeness” or “coverage” with native Julia methods. The quick summary of DifferentialEquations.jl 3.0 is the following: for hard problems, we now have methods specifically suited for the problem. We have methods for stiff ODEs, SDEs, DDEs, etc. and these work with the differential-algebraic forms of each of these equations. We still need to round out the suite, but I am pleased to say that for hard problems which require special methods that can be difficult to implement, we do have options available for you. Let’s go into some details.

Solvers for Stiff ODEs and DAEs

This is probably the area that will impact the most individuals. In DifferentialEquations.jl 3.0 we are happy to announce the release of a vast array of methods for solving stiff differential equations. While before we had wrappers for methods like CVODE from Sundials, LSODA, and radau from Hairer’s software, our offering here wasn’t too unique. However, we now have a wide array of high order methods for solving stiff ODEs. The centerpiece here are Rosenbrock methods and (E)SDIRK methods.

Rosenbrock methods are methods which are generally very good at lower tolerances. Hairer’s second book showed that high order Rosenbrock methods tend to be the most efficient methods when the required error is less than something around 5 digits. This is huge because this is the amount of accuracy many people want. Our new offerings of ODE solvers includes pretty much every Rosenbrock method that we could find that has been proposed in the literature. These methods have special interpolations so that sol(t) not only acts as a continuous function of the solution, but this continuous function is in some sense “stiffness-aware” and can with high accuracy produce the solution and its sharp turns between the solver’s steps. Being a “generic” Julia implementation, these all work with a wide array of Julia-defined number types, including high-precision arithmetic and (if the Jacobian is defined, see below) complex numbers. As far as we know, this is the first set of stiff ODE solvers with this flexibility. They all fully conform to the framework of DifferentialEquations.jl, meaning that they have event handling, the integrator interface, and all of the other extra goodies. They have adaptive timestepping with automatic initial dt calculations and all of the other features to make them a “fully automatic” solver. Using the mass matrices, these solvers can also handle DAEs.

In addition to the Rosenbrock methods we released a large set of (E)SDIRK methods. These methods have a specialized quasi-Newton solver for their implicit equations, making them highly efficient, especially in the case where there are times in which the Jacobian changes less quickly (it will skip factorizations when it determines that it can). As with the Rosenbrock methods, these are “fully automatic” and work with all of the event handling, generic numbers, etc.

So, how did we do? What do the benchmarks look like? DiffEqBenchmarks.jl is how we’ve been tracking some of our progress. In the benchmarks there we show that in the common test problems for stiff ODEs, these newest methods are the fastest we have available, even faster than CVODE from Sundials and radau from Hairer, in the “range of reasonable tolerances”, i.e. where the user wants an the error to be in the 9th digit or lower. I think that satisfies most uses cases and so we are pretty happy with the results. There are many other tests from users which report similar results that our new Rosenbrock methods and (E)SDIRK methods benchmark as the fastest for achieving the desired accuracy.

We are not surprised though: multistep methods like CVODE are specialized to decrease the number of function (f) calls which in turn is only useful when the system is sufficiently large. And fully implicit methods make use of larger linear systems to be efficient when more steps are required. Thus in the case where the system function is very costly or the number of ODEs is huge, Sundials will still be a good choice. Or if you need really high accuracy, radau will still be a good choice (note: these methods are still wrapped so you can keep using them). But other than these cases, we find our new methods to perform really well.

Let me mention a few caveats here which are left. One of them is complex numbers. Complex number handling is a little bit spotty in Julia packages right now. DifferentialEquations.jl fully supports them in each of the *DiffEq (OrdinaryDiffEq.jl, StochasticDiffEq.jl, DelayDiffEq.jl) solvers, along with all of the addons. However, where we run into issues is when interfacing with other packages. For example, ForwardDiff.jl and Calculus.jl cannot handle complex numbers. This becomes an issue only for the stiff solvers because the stiff solvers require the ability to calculate Jocobians, for which we use these packages. This is why I added the caveat “if you provide your own Jacobian”. But actually, we get the capability to numerically compute Jacobians with complex numbers by replacing Calculus.jl with DiffEqDiffTools.jl in DifferentialEquations.jl 3.0 (a little pre-mature, but the PR is just about done here). But in a week or so this will be a problem of the past, and we hope to integrate our tools into things like Optim.jl and NLsolve.jl so that more packages support complex numbers (you quantum physicists keep emailing me! ).

Solvers for Stiff Delay Differential Equations

Our new methods in OrdinaryDiffEq.jl, the high-order Rosenbrock and (E)SDIRK methods, extend over to stiff delay differential equations. These are specialized so that way they can “re-use” step information to be more efficient than classic designs. We only know of one other available free stiff delay differential equation solver (Radar5), but since we couldn’t find out how to get it to work (it requires some really intricate compilation binding so I don’t think it can be wrapped) we don’t have anything to benchmark against. But from what we’ve seen it works well! Once again (as always), DiffEqBenchmarks.jl is the open resource for seeing how things are doing.

Solvers for Stiff Stochastic Differential Equations

We wouldn’t be complete without saying that we also have methods for stiff stochastic differential equations. These are based on the SDIRK architecture of the ODE solvers and thus employ the same tricks to get efficiency. Not much more to say here.

Solvers for Ordinary, Stochastic, and Delay Differential-Algebraic Equations

Many of our stiff solvers allow for defining a mass matrix. The mass matrix is allowed to be singular, in which case the stiff solver will solve a differential-algebraic version of the ODE, SDE, or DDE. As far as we know, this is the first available set of solvers for SDAEs and DDAEs.

Solvers for for Second Order ODEs

Okay, you can convert first-order ODEs to second order ODEs and solve them like that. However, when doing so you don’t make use of the full structure of the second order equation, and thus you don’t get the full efficiency out. Runge-Kutta Nystrom methods are made directly for second order ODEs and we now have these methods implemented. In addition, in many cases one may want to solve an equation in a way that you know certain quantities are preserved over long-time integrations. These methods are known as symplectic integrators, and we have implemented a large array of symplectic integrators.

The format for these is what we call a “dynamical ODE”. The basic way to specify a DynamicalODEProblem is by specifying it as a second order ODE. However, we also allow one to directly specify the Hamiltonian for a physical system from which we use autodifferentiation to derive the equations of motion. Additionally, we allow a partitioned ODE form which allows one to specify the velocity component directly and thus allows for more advanced dynamics than a simple second order ODE. All of these problems can also directly be solved by first-order ODE solvers which will automatically do the conversion.

Solvers for State-Dependent Delay Differential Equations

State-dependent delay differential equations are delay differential equations where the delay factor depends on the differential equation itself. For example, you can say something like, the amount of growth in the fish population is dependent on the population from a few weeks ago (since that’s when conception would have occurred), but when there’s more fish there’s a longer delay since the development process is slowed when resources are scarce. This means that the derivative of now is dependent on the derivative of the past, but how long in the past is dependent on the value right now!

Delay differential equations are complicated to solve because these delays propagate discontinuities. If you don’t properly handle the discontinuities then you will not achieve high accuracy. For constant-time delays you can know exactly where all of the discontinuities will be a priori, and thus you can have the solver hit exactly those points in time in order to not have issues. For state-dependent delays, the timepoints for the discontinuities are dependent on the solution itself, so you need the numerical solution in order to know how to handle the discontinuities!

If you ever step over a discontinuity, you will suffer from increased error. Or do you? Questioning this assumption gives you the residual control methods. These are adaptive methods which have a robust form of error estimation and thus try to detect discontinuities by stepping over them and seeing if the resulting error is high. This is the method that MATLAB’s ddesd uses, and thus we implemented this as well. However, shortly before doing so, I received an email from some numerical delay differential equation researchers who questioned the validity of this approach because they ran MATLAB’s ddesd on some test problems and found out that its error was quite high. Well, our resdiual control methods match this behavior: they don’t tend to get more than 3 digits of accuracy, but they are pretty fast. To be fair, Shampine’s paper on ddesd said that it was for getting plotting accuracy, and not necessarily scientific computing accuracy.

So that method handles one case, what about the high accuracy case? The JuliaDiffEq contributor David Widmann is who to thank for this. Using the event handling setup in the ODE solvers, we setup a system by which the solver would continuously track and detect discontinuities, and use this to pullback and hit discontinuities “exactly”. Testing against numerical solutions this method is able to get to full floating point accuracy. This method is also compatible with all ODE solvers via method of steps, and thus allows for using stiff solvers and differential-algebraic delay equations via mass matrices.

Solvers for Boundary Value Problem (BVP) Solvers

This is a result from which you can thank Google Summer of Code. Boundary-value problems are extensions of ODEs which allow you to set conditions which the ODE must satisfy. Normally one thinks of the two-point boundary value problem where these conditions specify values that the ODE must be at the start and the end of the solution interval. We did create a method for two-point BVPs which mirrors that of bvp4c (though adaptivity is coming soon), but we generalized the allowed BVPs quite a bit. For many of the methods, you are able to specify conditions using the full solution and its interpolation. Thus one can make a “boundary condition” that the maximum of the second derivative over the full interval is 1. We honestly do not know of problems which utilize this full generality yet (though I do know of “multipoint BVP problems” which are of course a subset of this), so we’d like to hear if you end up using this for something crazy.

Partial Differential Equation (PDE) Toolkits: Linear Operators and FEniCS

Oh PDEs. If you watch my JuliaCon workshop, you’ll see that the same two questions always come up: what about stiff solvers, and what about PDEs? I just told you about solvers for stiff differential equations for a wide variety of problems, and now lets address our PDE tools.

Early on in DifferentialEquations.jl I created a finite element toolbox to go along with the software. It was very basic, and I realized that approach will not scale, so instead what we decided to do was wrap the popular FEM library FEniCS. This was part of a Google Summer of Code project which created FEniCS.jl. You can read the blog post which introduces it and what’s cool is that the pieces that FEniCS created, the assembled operator equations, can be directly converted to sparse matrices in Julia which can be used to solve time-dependent PDEs in our ODE solvers (or time-independent problems can just solve the implicit equation using whatever linear solver you choose from Julia).

But not every problem needs finite element methods. To help with finite difference methods, another GSoC project developed DiffEqOperators.jl which makes it easy to discretize PDEs via finite difference methods. Essentially, you tell it the derivatives you want to discretize and it spits out lazy (matrix-free) linear operators which are fully multithreaded and perform the stencil. Once again, this makes it easy to define the discretized ODE system from the PDE and then solve it using the ODE solvers. We also include upwind operators for stable discretizations of hyperbolic PDEs.

As you can see, this is a toolbox for solving PDEs. People who have a little bit of prior knowledge in solving PDEs can easily use these tools to build a method that solves their specific PDE. However, what we plan to do next is to use this toolbox to make some pre-made solvers for some common PDEs like diffusion-advection equations. With this development it will really complete our PDE story.

Conclusion: DifferentialEquations.jl 3.0 addresses the major concerns of the past

The main conclusion is this: people wanted methods for all sorts of solvers for stiff ordinary, stochastic, and delay differential equations, along with the differential-algebraic and partial differential variants, and they wanted these to work with generic numbers, event handling, all of the addons, etc. This announcement’s tl;dr is simply that we listened, and we released. Of course, we aren’t done (there’s always more to do), but what we can say is that it is highly likely that one of our offering will solve your differential equation well.

So what’s next? Well, we can always add some more methods which handle specific special cases better. That’s the goal of DifferentialEquations.jl 4.0 (and beyond). Here’s a look into what we have planned.

Multistep Methods (Adams, BDF) and Implicit ODEs

Classic multistep method solvers like LSODE and CVODE are some of the most commonly used methods. In most cases they aren’t the most efficient: this is a fact noted in Hairer’s benchmarks and now in ours. However they have major upsides when the user’s function f is expensive, or when the system of ODEs is large. This is something that comes up in large PDE discretizations and is why these methods are central to solving large-scale PDEs. We have put this off because it is a niche area and we have pretty good wrappers to the classics like Sundials, LSODA, DASKR… but it is definitely time that we tackle these methods with a native Julia implementation.

One other thing to mention here is that multistep BDF methods are also what have traditionally gave rise to fully implicit DAE solvers like DASSL, DASPK, and IDA. This is an area where we have been lacking quite a bit in terms of native solver capabilities (mostly relying on wrappers), so we will need to spend some time here. We also need to build tooling for finding consistent initial conditions like is done in these solvers… we’re a ways off here. Modeling tools like Sims.jl and Modia.jl directly utilize IDA since we don’t offer anything of interest in this area, but hopefully we can develop some native Julia tooling here and help link it to these other packages so that we can expand the capabilities of not only us but modeling packages as well. If we can provide a better solving backend and they provide a great modeling front end, Julia will be the star of this field.

Fully implicit Runge-Kutta methods also have a niche. Some methods, like radau, are great for high-accuracy (low tolerance) solving of stiff equations. Others are high order symplectic methods for stiff differential equations. They also do really well with DAEs in mass matrix form. Now that we have the availability, these are definitely areas that we will tackle in the near future as well, not just in ODEs, but also in SDEs (and note that the ODE part builds DDE solvers for free as well).

Exponential Integrators

Exponential integrators allow you to exploit linearity in the definition of an ODE, SDE, or DDE. There are two forms of interest: $u'=A(t)u$ for $A$ is a time-dependent linear operator, and $u'=Au + f(t,u)$ where $A$ is a time-indpendent linear operator. The first form shows up in a lot of quantum mechanics situations. The latter comes from discretizations of semilinear PDEs. Both of these can be solved with standard first-order ODE solvers, but the efficiency can be improved by using $A$ directly.

We have already made great strides in this direction. There are some solvers released for both types of equations, and we have developed an interface, the DiffEqOperator, for handling the definition of $A$ in a way the solvers can exploit. However, the crucial linear algebra tools were picked up by Marcelo Forets and implemented in ExpoKit.jl, and using their expmv and phimv implementations we can tackle the higher order methods. I wouldn’t expect this until the next summer since I see portions of this project as a great Google Summer of Code project, so if you’re interested please feel free to get in contact with us.

Implicit-Explicit (IMEX) Methods

IMEX methods, where the user can split the function f into two portions so that way one part is explicit and one part is implicitly solved (i.e. a stiff and nonstiff part) has recently received lots of popularity for solving PDEs. We have most of the pieces for high order IMEX methods. The ESDIRK methods from Kennedy and Carpenter are the additive Runge-Kutta methods that are used in Sundials’ ARKODE IMEX solvers. We just haven’t added the explicit part.

But we have the architecture which allows the user to define IMEX methods. There’s also plenty of other IMEX methods which can be implemented as well. Since the architecture here is “already done”, it’s simply a matter of coding in a the inner loop for a few new methods. To me, this sounds like a great Google Summer of Code project as well, so we may be holding off on development here until the next summer.

Finding Out Our Partial Differential Equation (PDE) Interface

I was happy to release (at least the beginnings of) PDE toolkits… but that satisfies a small group of people who know numerical methods for PDEs and want these pieces in order to write solvers more easily. In practice, many scientists probably don’t know how to (or don’t want to) do this (they are specialists in science! Not solving PDEs!). We need to provide something like MATLAB’s pdepe: simple interfaces for solving common PDEs. While one way we will build these solvers will be to use our toolbox and build method of lines integrators, we will need to make use of the distributed architecture of DifferentialEquations.jl in order to get good coverage of the PDE landscape. Indeed, I know of some individuals like John Gibson who are building spectral PDE solvers and really seem to know what they are doing, and so it would be a shame if we didn’t have a way to allow users to directly interface with these tools (it’ll also be a great way for methods researchers to easily benchmark, hopefully pulling an even larger community of developers in).

The answer will be some kind of problem hierarchy where the user defines something like a DiffusionAdvectionProblem where it contains a bunch of functions and there is a common interpretation of what the problem definition is and how solvers should treat it, and then they should all use that to spit out a similar solution. We have the pieces of how that can work in our (fantastically outdated) FEM Heat and Poisson methods, but the issues are finding out how to do things like boundary conditions in a way that we can make the most out of every’s PDE solver while not complicating the common interface. Once we have solvers all wrapped together, there will still be a nightmare: how do you document this? If we have a different set of four packages available for 5 different PDEs, the combinatoric explosion in documenting the nuances means we can’t add all of this to our current documentation (which is already huge). So we may need a “different section” of the docs somehow? To me it’s very unclear how we will document this well, so my plan is to just start adding the functionality and find out how to document it as we go along. It will probably need to be re-written a few times before it’s any good, so please bear with us!

High Order Adaptive Solvers for Stiff SDEs

This is actually one of my recent research projects. I have new methods for high order adaptive solvers for stiff stochastic differential equations, and will be submitting the publication soon. When this is published, the associated methods will be released and we will have high order adaptive methods for stiff SDEs. So stay tuned!

Flesh out the BVP solvers

Our Shooting methods are very flexible, but they will never do well on problems which are sensitive to the initial condition. We need to instead make our MIRK methods get all of the bells and whistles: continuous extension, adaptivity, etc., and we need to wrap some of the classic Netlib solvers into the same interface so we can start to do some comprehensive benchmarking.

Parallel-in-time ODE Solvers

During the last Google Summer of Code we took a stab at developing some parallel-in-time ODE solvers using Neural Networks. While the student did a great job at trying a bunch of different strategies, we have come to realize that neural networks simply do not “know” enough about the structure of the differential equations in order to be efficient. When talking to a friend about PDE solvers, he noted that he tested the efficiency of TensorFlow’s neural net PDE solver (which it has in a tutorial) against more standard methods and noted similar issues with efficiency. So we tried, and it was definitely a very interesting research project, but this direction didn’t yield the results that we hoped.

Thus instead I am hoping to take a step over to different methods. Parallel-in-time integration methods like parareal integration and the XBraid software do exist, and so I plan on taking a stab on adding these to the repertoire of DiffEq. Note that for large and expensive ODEs, SDEs, and DDEs, one can already parallelize the calculation of f using the current tooling. Thus these methods are for when you want to solve problems where the set of ODEs is quite small, yet you need to solve over a large timespan and have a lot of parallel computing power available. So once again, there’s no rush as this is quite a niche, but we plan to get to it in the next release cycle. In fact, I hear it’s so niche that I was told in an email that parareal is only good for long time integration when you have >128 cores available… let’s make an open-source implementation and try it out ourselves.

Automatic Stiffness Detection and Switching

What’s a stiff differential equation? That can be hard to explain and predict. For this reason, many people like the idea of automatic stiffness detection and allowing the method to automatically switch between solvers to handle the different types of equations more effectively without user input. We have all of the tooling for building this, and the user can actually specify switching strategies themselves using our CompositeAlgorithm setup, but in the next release cycle I hope to release some methods which have this built in, akin to something like LSODA. Note that my newest paper on SDE methods includes stiffness detection as well, so in one fell swoop we plan to add stiffness detection and switching for ODEs, SDEs, and DDEs (and of course differential-algebraic version via mass matrices, and …, you get the picture about how all of this stuff is composed together!)

This is actually at a point where it would not be too difficult to do but would take a good chunk of time and also take some research time, so I am putting it off and hoping this can be a Google Summer of Code project or another project with a student.

Efficiency Improvements for Cheap DEs

One final wart in the DiffEq architecture. Due to issues with Julia, it seems that if an ODE is “sufficiency cheap”, i.e. takes less than a millisecond to solve, our setup has some inefficiencies. This is from how we do the setup of the integrator and limitations Julia has on type inference. This issue is constant, meaning that for more expensive ODEs/SDEs/DDEs it’s still just a millisecond. However, it’s annoying and we hope to remedy this. In reality, we haven’t seen this affect real-world benchmarks since < millisecond ODEs are not necessarily where people are looking at performance, but it does start to become an issue when attempting to do something like parameter estimation on cheap ODEs since this involves solving the same ODE thousands of times. Part of the issue is due to the use of keyword arguments whose performance will be fixed in Julia v0.7. For the inference issues, we hope to get the right updates into an early Julia v1.x, along with making a few structural changes, and then this issue will go away. Other enhancements along this line will be a reinit interface to make it easier to reuse internal caches and cleaning up the solution type pre-allocation interface.

Conclusion

At this point, we are very satisfied with our offering. The 30 people who make up the JuliaDiffEq team have really built a software which has the methods to solve “most” differential equations that users encounter, and also do so efficiently. In the coming months we hope to add extra methods for specific and important niches to our offerings and fill in some holes. But together, I think we have a pretty solid offering, and everything else is (important) icing on the cake.

The post DifferentialEquations.jl 3.0 and a Roadmap for 4.0 appeared first on Stochastic Lifestyle.