In this and the next chapter on Useful Packages I want to give an overview of the functionality and scope Julia can deliver. Details you can find in the referenced links.
# Packages we need here:
import Pkg; Pkg.add(["PyCall"])
There are plenty of Julia datatypes I have not yet mentioned. This includes dictionaries and tuples:
# Construct a dictionary, note that dictionaries are typed as well.
a = Dict("a" => 1, "b" => 2)
println(typeof(a))
for (l, i) in a # Loops are over key-value pairs
println("key $l maps to $i")
end
# Does not work:
a[2] = 5
# A tuple
t = (2, "abc", )
@show t
# A named tuple
nt = (value=2, description="abc")
@show nt
println("The description is ", nt.description)
Apart from the built-in types, Julia also offers to declare custom types. Both abstract types and concrete types can be defined:
abstract type MyAbstract end
struct MyPoint <: MyAbstract
x::Float64
y::Float64
description # Note: Fields do not need to be typed, in this case they are Any
end
"Class methods" or "member functions" as such do not exist in Julia. Instead one uses free functions with type-specific dispatch, e.g.
import Base: length
length(p::MyPoint) = sqrt(p.x^2 + p.y^2)
shift_in_x(p::MyPoint, v) = MyPoint(p.x + v, p.y, p.description)
p = MyPoint(3, 4, "bla")
length(p)
Compared to other languages such structs are used a lot less and much functionality in the standard library and third party packages avoids them completely. Still they can be very useful to combine multiple enties, which belong together.
Julia comes with a built-in library for file system interaction. This makes a lot of use for the do
syntax we saw before. For example:
cd(homedir()) do
println("All files and directories in your home are:")
println(readdir())
end
# Back to the old directory
# Open a file for writing:
open("testfile", "w") do f
write(f, "test")
end
# Open a file for reading:
open("testfile") do f
println(readline(f))
end
rm("testfile")
Some functionality for unit testing is part of Julia's standard library. As of now it is pretty basic, but it get's the job done:
using Test
@testset "My first test" begin
@test 1 == 1 # Should be good, no?
@test_throws BoundsError (1, 2, 3)[5]
end
Julia has pretty strong metaprogramming support. In fact macros are even an integral part of the core language and their use is frequent in the standard library itself. Common use cases include a temporary change of behavior, to automate the generation of repetitive code or to hide some implementation details under a more convenient syntax. In this course we have already met quite a few macros, for example @show
, @time
, @btime
, @.
, @gif
, .... It is extremely uncommon that one needs to add or modify macros for daily work with Julia. Still if you ever come across the need to implement a macro, their syntax is quite readable:
macro make_addfunction(value)
:(# Start an expression, which is evaluated as code
# $(Symbol(:add, value)) makes a code symbol by concatenating
# two others together
function $(Symbol(:add, value))(x)
x + $value
end
)
end
@make_addfunction 10
add10(4.5)
As discussed before type-stability is important in performance-critical parts of the code such that Julia can infer concrete types before JIT compilation and be very specific on the emitted machine code.
But introducing type instabilities can come easier as you might think. Consider the following unsuspicious piece of code for example:
# Compute the sum of a column
function sumcol(a, n)
@assert ndims(a) == 2
accu = 0
# @inbounds disables bounds check for a[i, n]
@inbounds for i in 1:size(a, 1)
accu += a[i, n]
end
accu
end
function average_sumcol(a)
@assert ndims(a) == 2
res = 0
for j in 1:size(a, 2)
res += sumcol(a, j)
end
res / size(a, 2)
end
using BenchmarkTools
a = randn(ComplexF64, 500, 500)
@btime average_sumcol(a);
It's reasonably fast, but it actually contains subtle type instabilities. This can be found using @code_warntype
:
@code_warntype average_sumcol(a);
This is a rather complicated output, but it roughly speaking tells us with the bold print, that somehow the compiler cannot decide between Int64
and Complex{Float64}
for some lines of code and not between Float64
and Complex{Float64}
in others. The variables blocks tells us that the res
variable declared in line 14 might be a problem. And indeed the issue is that we declare it as = 0
, which makes it an integer. Therefore the return type of average_sumcol
is not stable, because it is different for empty arrays (Int64
) and non-empty arrays (Complex{Float64}
). A similar spot is line 4. Replacing both instances with zero(eltype(a))
makes the code completely type stable (and actually a little faster).
In this example one might still be able to read and understand the output, but especially for larger functions @code_warntype
can become long and unclear. Of help is a package with the unrememberable name Cthulhu
and its macro @descend_code_warntype
. Unfortunately it does not work in Jupyter notebooks at the moment, so if you want to see it, try it in the REPL.
In this section, we want to very briefly explore one (of many) ways to do multithreading in Julia. For this we first need to increase the number of threads the Julia kernel uses. We check on the CPUs of this machine
Sys.cpu_summary()
In my case, I can set the number of threads to 4
:
Threads.nthreads() = 4
Let us recall the computation of the total potential for the multidimensional case in Dancing_Particles:
using LinearAlgebra
Vmorse(r; re=0.7, α=1.3, D=10) = D*(1 - exp(-α * (r - re)))^2 - D
function Vtot(Vpair, x)
n_particles = size(x, 2)
accu = eltype(x)(0) # Get a zero of the appropriate type
for i in 1:n_particles-1, j in i+1:n_particles
accu += @views Vpair(norm(x[:, i] .- x[:, j]))
end
accu
end
We can test it on 5000 particles:
using BenchmarkTools
x = randn(2, 5000)
@time Vtot(Vmorse, x)
Puh ... that's getting slow. An obvious place to parallelise is the accumulation over particles. For this we will use the Threads.@threads
macro. The other required change is to use an atomic accumulator to avoid inconsistencies:
import Base.Threads: @threads, Atomic
function pVtot(Vpair, x)
n_particles = size(x, 2)
T = eltype(x)
accu = Atomic{T}(zero(T)) # Atomic accumulator
@threads for i in 1:n_particles-1
local_accu = zero(T) # Thread-local accumulator
for j in i+1:n_particles
local_accu += @views Vpair(norm(x[:, i] .- x[:, j]))
end
accu[] += local_accu
end
accu[]
end
@time pVtot(Vmorse, x)
Of course this is still far from perfect, since the load between the threads is not exactly balanced, but we'll leave it here.
For distributed computing over mulitple nodes, Julia offers DistributedArrays
, a package for spreading array data across workers and Distributed
from the standard library for classic building blocks of distributed computing (parallel maps, reductions, channels ...)
Mixing Julia and C, C++ or Fortran works very well. The underlying reason is because Julia gets compiled via the LLVM compiler infrastructure, which is capable of the other languages out of the box.
For example we can directly call the dot product function ddot
from LAPACK:
# Call to lapack to compute dot product
v = [1.0, 2.0, 3.0, 4.0]
w = [2.0, 3.0, 4.0, 5.0]
VdotW = ccall((:ddot_, "liblapack"), Float64,
(Ref{Cint}, Ptr{Cdouble}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}),
length(v), v, 1, w, 1)
println("v ⋅ w = $VdotW")
(actually ... this should not be done like this in practice, because even if one wants to avoid the builtin julia function dot
, there are much better LAPACK wrappers in LinearAlgebra.LAPACK
).
As another example, take the following piece of C code:
code = """
double sum_array(double* array, int n) {
double accu = 0.0;
for (int i = 0; i < n; ++i) {
accu += array[i];
}
return accu;
}
""";
Which we first compile using cc
into the shared library file libsums.so
:
open("sums.c", "w") do f
write(f, code)
end
run(`cc -shared -o libsums.so sums.c`)
rm("sums.c")
Using sum_array
from Julia now similarly only takes an appropriate ccall
:
v = [1.0, 2.0, 3.0]
res = ccall((:sum_array, "libsums"), Cdouble,
(Ptr{Cdouble}, Cint), v, length(v))
Of course compiling each time is not so convenient. If you require mixing C++ and Julia a lot, have a look at the Cxx
package. As a bonus it even adds an interactive C++ shell to Julia's REPL ...
By far the most advanced foreign-language integration has been achieved with Python. Here you can easily use any python-package from Julia. For example NumPy:
using PyCall
np = pyimport("numpy")
np.random.rand(3, 4)
A lot of type conversion happens in the background automatically. In this case NumPy arrays are automatically converted to Julia (and thus from row-major to column-major ...). The reverse (i.e. using Julia from python) is also possible, see Using_Julia_From_Python for a quick demonstration.
... you can find in the Julia documentation ;)