The Julia Language

08/17/19

Numerical scientific computing with Julia

Throughout my career I have learnt (to various degrees) the computer languages of (in chronological order): Matlab, C Shell, Bash, IDL, Fortran, C, Python, HTML, R, NCL, Grads, Hadoop, d3, SAS, Excel, SQL, Scala and Apache Spark. It is a natural process that new computer languages emerge as technology changes, attitudes change and funding changes. For example, I was taught Matlab during my undergraduate. During my PhD we had an IDL guru who would help the grad students. When I moved to Miami for my post-doc I had a clean slate. I opted to go for python as my research group used and it was becoming the emerging language of ocean and atmospheric science. Since then (2015) python has exploded (Figure 1). It is open source (free!) and it has a great community of people who are willing to go out of their way to help. For example, Matt Rocklin (NVIDA) was a great mentor for me working on dask-jobqueue.

Figure 1: Percentage of questions asked on StackOverflow. Taken from https://insights.stackoverflow.com/trends?tags=r%2Cpython%2Cmatlab

However, I like to think of my experience with other languages as other tools in my toolkit. Matlab has some robust probability distributions functions which was helpful when I worked at BP on offshore platform design criteria [1]. Fortran is still the gold standard when it comes to computation speed for scientific computing because of its ease of parallelization. I'm starting to enjoy the scalability of Apache Spark due its integration with the cloud computing platform databricks.

It is only a matter of time before another contender enters the ring. As much as this frustrates people: "i've just learned one language..." I see it as a good thing as it offers another tool for my toolbox. Other languages can also benefit from the new language. An example of this is Matlab taking inspiration from python for more aggressive broadcasting for operands of binary operators [2].

Enter the Julia language: Created in 2009 with the aspiration to be high-level and fast. Here is a good short blog by the creators on why they created Julia:

"We want a language that’s open source, with a liberal license. We want the speed of C with the dynamism of Ruby. We want a language that’s homoiconic, with true macros like Lisp, but with obvious, familiar mathematical notation like Matlab. We want something as usable for general programming as Python, as easy for statistics as R, as natural for string processing as Perl, as powerful for linear algebra as Matlab, as good at gluing programs together as the shell. Something that is dirt simple to learn, yet keeps the most serious hackers happy. We want it interactive and we want it compiled."

Lofty goal. Years down the line Julia is now at version 1.1 with 4 million downloads and 2,400 community developed packages. It is time I start catching up with the times. My interest is in numerical scientific computing and a couple of packages have recently caught my attention. One is on differential equations (which are at the heart of machine learning algorithms) and the other on time series prediction.

Installing Julia

These are instructions for my PC which runs Ubuntu:

$ wget https://julialang-s3.julialang.org/bin/linux/x64/1.1/julia-1.1.1-linux-x86_64.tar.gz
$ tar xvzf julia-1.1.1-linux-x86_64.tar.gz
$ sudo ln -s /home/ray/src/julia-1.1.1/bin/julia /usr/local/bin/julia

Installing IJulia

$ julia
julia> using Pkg
julia> Pkg.add("IJulia")
julia> exit()
$ jupyter notebook

Julia 101

For getting started with Julia there is https://julialang.org/learning/. I'm going to use the tutorials in here. Below are some of my notes; more comprehensive notes are kept here (including bench-marking):

  • Gets docs of print by doing ?println
  • Run shell commands by doing ;ls
  • println("Hello world!")
  • Julia figures out datatypes typeof(42)
  • You can program with emojis... \:smiley_cat: then <tab> then <enter>
  • fibonacci = [1, 1, 2, 3, 5, 8, 13]; push!(fibonacci, 21); pop!(fibonacci)
function sayhi(name)
    println("Hi $name, it's great to see you!")
end


function f(x)
    x^2
end
map(f, [1, 2, 3])

map(x -> x^3, [1, 2, 3])


broadcast(f, [1, 2, 3]); f.([1, 2, 3])
  • Create a 3x3 matrix of 1:9 A = [i + 3*j for j in 0:2, i in 1:3]; f(A) is dot product (A .A); f.(A) squares the elements
Pkg.add("Colors")
using Colors
palette = distinguishable_colors(100)

DifferentialEquations.jl

Taken from the docs: This is a suite for numerically solving differential equations in Julia. This packages solves Ordinary differential equations (ODEs), Partial differential equations (PDEs) and more. Here is an example of solving the equation du/dt = alpha * u. The analytical solution is u(t) = u0 * exp(alpha * t). Note if alpha is a negative number this gives an equation known as Newtonian cooling which is the rate your coffee cools at!

using Pkg
Pkg.add("DifferentialEquations")
using DifferentialEquations

f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0) # period to solve over
prob = ODEProblem(f,u0,tspan)

sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

using Plots
plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")

Solving the Lorenz 1963 equations. Note: this equations gave rise to the 'Lorenz attractor'. It is also a good analogy for weather forecasting in that the physical system can be completely deterministic and yet still be inherently unpredictable based on imperfect knowledge about initial conditions.

function lorenz(du,u,p,t)
 du[1] = 10.0 * (u[2] - u[1])
 du[2] = u[1] * (28.0 - u[3]) - u[2]
 du[3] = u[1] * u[2] - (8/3) * u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob)

plot(sol,vars=(1,2,3))

TimeseriesPrediction.jl

This package models time series using methods of nonlinear dynamics. In the example below it predicts a time series using a local method with two parameters (gamma and tau) and using 3 nearest neighbor points. eq

using Pkg
Pkg.add("TimeseriesPrediction")
Pkg.add("DynamicalSystemsBase")
Pkg.add("PyPlot")

using TimeseriesPrediction # re-exports DelayEmbeddings
using DynamicalSystemsBase # to access some systems

ds = Systems.roessler(0.1ones(3))
dt = 0.1
data = trajectory(ds, 1000; dt=dt)
N_train = 6001
s_train = data[1:N_train, 1]
s_test  = data[N_train:end,1]

ntype = FixedMassNeighborhood(3)

p = 500
s_pred = localmodel_tsp(s_train, 4, 15, p; ntype=ntype)

using PyPlot
figure()
plot(550:dt:600, s_train[5501:end], label = "training (trunc.)", color = "C1")
plot(600:dt:(600+p*dt), s_test[1:p+1], color = "C3", label = "actual signal")
plot(600:dt:(600+p*dt), s_pred, color = "C0", ls="--", label="predicted")

title("Pool of points: $(N_train), predicted points: $(p)")
xlabel("\$t\$"); ylabel("\$x\$")
legend(loc="upper left")
tight_layout()