Dynamical systems in Julia
Introducing Julia
Julia is a powerful language for scientific computing. It is a high-level language that is a lot about performance. Its main feature is the use of a just-in-time compiler to achieve a speed for numerical calculations that is close to C or Fortran code. The learning curve can be a bit steep, especially since one needs to have a good understanding of the type system to work efficiently in Julia. That said, it was something like love on first sight.
These past days I have been implementing a simple ice structure interaction model, that is basically a system of ODEs with discrete events that change the system state randomly. This is a non-smooth dynamical system, a very interesting class with many fascinating properties. I spare you the details here!
Implementing the model and running it was quite straightforward, but I had a number of challenges with Julia itself and the bifurcation analysis. The latter relies on PyDSTool, a Python library that is called from within Julia.
My first tests were with JuliaPro, a version of Julia that comes with an IDE and a comprehensive set of pre-installed libraries, similar to Anaconda for Python. This is a nice idea for getting productive quickly with Julia, but for this project it turned into a nightmare. For some reason the IDE (Juno, building on the Atom editor) was using a different directory for the compiled Julia code than the command-line version of Julia, and it was constantly recompiling my functions. Something was definitely not working right, and I had a lot of trouble getting the DifferentialEquations.jl package compiled. Some parts of it did not compile, and I figured out that it was missing dependencies. But some other issues were a complete myster to me. In the end I gave up, removed JuliaPro and tried to install everything myself. This was much easier than I thought, and worked out beautifully in the end.
Here I want to share my setup, it might be useful for someone else out there. And besides, I keep referring back to these notes when I need to setup Julia on another machine.
Scope
The following assumes a 64bit Linux system. I have tested the setup on two different machines with Linux Mint, so it should work on e.g. Ubuntu with no modifications. I will soon test it on a workstation with Sabayon Linux, too.
Let’s get going
The first thing we do is download the official Julia binary release from julialang.org. Currently this is the 0.6.1 version of Julia. It is unpacked into a subdirectory of our home, e.g. $HOME/local
I assume Python to be present in the system, in my case this is currently Python 2.7.12. Let us install ipython and jupyter, or upgrade them to the latest version, if already existing.
$ sudo -H pip2 install --upgrade ipython jupyter
Note that I am using the pip2 command here, to make sure that the packages are installed for Python 2.x, not for Python 3.x. It should be already installed with your Python distribution, otherwise the package manager of your distribution should be used.
Then setup environment variables with the paths to the Python interpreter and the jupyter command.
$ declare -x PYTHON=/usr/bin/python
$ declare -x JUPYTER=/usr/local/bin/jupyter
Let us test that jupyter is working. If properly configured, the following command should fire up the browser of your choice with jupyter notebook running:
$ jupyter notebook
Fascinating, but let’s kill this server now (CTRL-C) and continue. For plotting inside the Atom editor (see below), we will also need the matplotlib library:
sudo -H pip2 install --upgrade matplotlib
Make sure that Tkinter is installed. In Mint, this would be the python-tk package.
Now we setup a few environment variables to tell Julia where to find its code (this should be the directory containing the julia executable) and packages (this defaults to ~/.julia/v0.6). Of course we also want to make sure that julia lies on our PATH.
declare -x JULIA_HOME=$HOME/local/julia/bin
declare -x PATH=$JULIA_HOME:$PATH
declare -x JULIA_PKGDIR=$HOME/.julia/v0.6 # default
The time has come to install the PyCall and IJulia packages from inside Julia, and support for plotting.
julia> Pkg.add("PyCall")
julia> Pkg.add("IJulia")
julia> Pkg.add("PyPlot")
julia> Pkg.add("Plots")
Note that specifying PYTHON in the environment (see above) means that PyCall will use the already installed Python (/usr/bin/python in this case). Let us give PyCall a testrun.
julia> using PyCall
INFO: Precompiling module PyCall.
julia> @pyimport math
julia> math.sin(math.pi / 4) - sin(pi / 4)
0.0
Want an IDE? Let us install atom. Download the .deb file, double-click on it, and voila.
We open atom, go to Edit -> Preferences -> Install and install the uber-juno atom meta-package. This will download and install a lot of goodies, so be patient. But we still do not have a Julia console in atom, although a Julia menu now has appeared in atom. If we try Julia -> Open Console and type “2 + 2” in the Console window, a whole lot of things will start to happen… until finally Julia runs for the first time inside atom and we get the expected answer “4”… I was not too happy with atom before, but this is really neat.
And while we are already busy with this editor, let’s install the project-viewer package, too. We can then open a Project-Viewer pane, create a Julia group in it, and add our project paths conveniently.
Let us try some plotting in atom now. We switch to cycler in the Julia -> Settings -> Boot Mode options first. Then the following should produce two plots, the first one in the atom Plots pane, the second one in a separate window:
julia> using Plots
julia> x = linspace(0,2*pi,1000); y = sin.(3 * x + 4 * cos.(2 * x));
julia> plot(x, y)
julia> workspace()
julia> using PyPlot
julia> plot(x, y, color="red", linewidth=2.0, linestyle="--")
Analysis of dynamic systems
Back to business: I want to solve ordinary differential equation models with Julia. The package to install is DifferentialEquations.jl.
julia> Pkg.add("DifferentialEquations")
Let us test our installation. We try the ODE tutorial from the DifferentialEquations documentation - I will not repeat it here. This should runs all fine, including the plots.
But we also want to do bifurcation analysis, for which we need the PyDSTool functionality in Python. PyDSTool is notorious for being not so straightforward to setup. Let’s see.
I try to install this Python library directly
sudo -H pip install --upgrade PyDSTool
which currently installs PyDSTool version 0.90.2 (it will also upgrade our numpy and scipy libraries to the latest version). Using PyDSTool however, results in a first surprise, a Runtime error:
SciPy v0.5.1 or above is required This is strange, as my scipy version is currently 1.0.0. Looking into the main init.py file of PyDSTool (/usr/local/lib/python2.7/dist-packages/PyDSTool/init.py) shows a bug! The check for the version number in line 75 only looks at the middle digit of the version number: if vernums[1] < 5: raise RuntimeError(“SciPy v0.5.1 or above is required”) What is needed is the following if vernums[1] < 5 & vernums[0] < 1: …
After this patch, we are not done, however. Importing PyDSTool
>>> import PyDSTool as dst
results in a mysterious error message
‘module’ object has no attribute ‘sph_jn’ connected to the scipy library, apparently.
The solution is to downgrade scipy to, for example, version 0.16.0, since the definition of the Bessel functions spn_jn was deprecated in scipy==0.18.0 (as one can read in the scipy release notes).
$ sudo -H pip2 install --upgrade --force-reinstall scipy==0.16.0 PyDSTool
Let us now do some testing. First the plotting:
>>> import matplotlib.pyplot as plt
>>> plt.plot([1, 2, 3, 4])
>>> plt.show()
Then we can try one of the PyDSTool tutorials, e.g., the Calcium channel model. If this runs, congratulations!
Now on to Julia itself. We might need to rebuild PyCall and IJulia:
> Pkg.build("PyCall")
> Pkg.build("IJulia")
Let us test PyDSTool with the Calcium channel model again, but in Julia this time.
