Other features worth mentioning

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.

In [ ]:
# Packages we need here:
import Pkg; Pkg.add(["PyCall"])

Dictionaries and tuples

There are plenty of Julia datatypes I have not yet mentioned. This includes dictionaries and tuples:

In [ ]:
# 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
In [ ]:
# A tuple
t = (2, "abc", )
@show t

# A named tuple
nt = (value=2, description="abc")
@show nt

println("The description is ", nt.description)

Custom types and structs

Apart from the built-in types, Julia also offers to declare custom types. Both abstract types and concrete types can be defined:

In [ ]:
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.

In [ ]:
import Base: length
length(p::MyPoint) = sqrt(p.x^2 + p.y^2)
In [ ]:
shift_in_x(p::MyPoint, v) = MyPoint(p.x + v, p.y, p.description)
In [ ]:
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.

Interacting with the filesystem

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:

In [ ]:
cd(homedir()) do
    println("All files and directories in your home are:")
    println(readdir())
end
# Back to the old directory
In [ ]:
# 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")

Unit testing

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:

In [ ]:
using Test
@testset "My first test" begin
    @test 1 == 1  # Should be good, no?
    @test_throws BoundsError (1, 2, 3)[5]
end

Macros and metaprogramming

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:

In [ ]:
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
In [ ]:
@make_addfunction 10
In [ ]:
add10(4.5)

Checking up on the compiler (@code_warntype and friends)

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:

In [ ]:
# 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
In [ ]:
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:

In [ ]:
@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.

Multithreading

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

In [ ]:
Sys.cpu_summary()

In my case, I can set the number of threads to 4:

In [ ]:
Threads.nthreads() = 4

Let us recall the computation of the total potential for the multidimensional case in Dancing_Particles:

In [ ]:
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:

In [ ]:
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 [email protected] macro. The other required change is to use an atomic accumulator to avoid inconsistencies:

In [ ]:
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
In [ ]:
@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.

Distributed computing

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 Fortran, C, python, R, Java, ...

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:

In [ ]:
# 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:

In [ ]:
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:

In [ ]:
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:

In [ ]:
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:

In [ ]:
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.

Other supported languages include Java or R, which can be used to import functionality from the respective languages into Julia.

More details

And even more features ...

... you can find in the Julia documentation ;)