# Speed by specialisation

In the previous notebooks we discussed Julia's elaborate type system and how multiple dispatch allows functions to decide on the method to dispatch really from looking at all the arguments. Intuitively **more specialised code**, i.e. code that is allowed to exploit known structure, **is faster.** In this notebook we will explore how Julia is able to exploit this while retaining generic code ... thanks to multiple dispatch.

## The speed of Julia

Now let us first try to address the question *Is Julia fast?*

To some extend this is a bit of a mismatching question, since one is able to write slow code in any language ... so let's try do address something else instead: *Can Julia be fast?*

## A simple example: Sums

Let's compare for a simple example ...


```julia
function mysum(v)
    result = zero(eltype(v))
    for i in 1:length(v)
        result += v[i]
    end
    result
end
```


```julia
v = randn(10^7);   # Large vector of numbers
```


```julia
using BenchmarkTools
```


```julia
# Compile some C code and call it from julia ...
using Libdl
code = """
#include <stddef.h>
double c_sum(size_t n, double *v) {
    double accu = 0.0;
    for (size_t i = 0; i < n; ++i) {
        accu += v[i];
    }
    return accu;
}
"""

# Compile to a shared library (with fast maths and machine-specific)
const Clib = tempname()
open(`gcc -fPIC -O3 -march=native -xc -shared -ffast-math -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, code) 
end

# define a Julia function that calls the C function:
c_sum(v::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(v), v)
```


```julia
bench = @benchmark c_sum($v)
```


```julia
times = Dict()
times["C (naive)"] = minimum(bench.times) / 1e6
times
```

An explicitly and fully vectorised C++ code:


```julia
# Compile some C code and call it from julia ...
using Libdl
code = """
#include <stddef.h>
#include <immintrin.h>

double c_sum_vector(size_t n, double *v)
{
    size_t i;
    double result;
    double tmp[4] __attribute__ ((aligned(64)));

    __m256d sums1 = _mm256_setzero_pd();
    __m256d sums2 = _mm256_setzero_pd();
    for ( i = 0; i + 7 < n; i += 8 )
    {
        sums1 = _mm256_add_pd( sums1, _mm256_loadu_pd(v+i  ) );
        sums2 = _mm256_add_pd( sums2, _mm256_loadu_pd(v+i+4) );
    }
    _mm256_store_pd( tmp, _mm256_add_pd(sums1, sums2) );

    return tmp[0] + tmp[1] + tmp[2] + tmp[3]; 
}
"""

# Compile to a shared library (with fast maths and machine-specific)
const Clib = tempname()
open(`gcc -fPIC -O2 -march=native -xc -shared -ffast-math -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, code) 
end

# define a Julia function that calls the C function:
c_sum(v::Array{Float64}) = ccall(("c_sum_vector", Clib), Float64, (Csize_t, Ptr{Float64}), length(v), v)
```


```julia
bench = @benchmark c_sum($v)
```


```julia
times["C (vectorised)"] = minimum(bench.times) / 1e6
```

How does our version do?


```julia
bench = @benchmark mysum($v)
```


```julia
times["Julia (naive)"] = minimum(bench.times) / 1e6
```

But unlike C we have not yet tried all tricks!


```julia
function fastsum(v)
    result = zero(eltype(v))
    @simd for i in 1:length(v)    # @simd enforces vectorisation in the loop
        @inbounds result += v[i]  # @inbounds suppresses bound checks
    end
    result
end

# Still nicely readable code ...
```


```julia
bench = @benchmark fastsum($v)
```


```julia
times["Julia (simd)"] = minimum(bench.times) / 1e6
```

Compare with python ...


```julia
using PyCall
numpysum(v) = pyimport("numpy").sum(v)

bench = @benchmark numpysum($v)
```


```julia
times["Numpy"] = minimum(bench.times) / 1e6
```


```julia
py"""
def py_sum(A):
    s = 0.0
    for a in A:
        s += a
    return s
"""

pysum = py"py_sum"

bench = @benchmark pysum($v)
```


```julia
times["Python (naive)"] = minimum(bench.times) / 1e6
```

Overview ...


```julia
times
```

## A more complicated example: Vandermonde matrices
(modified from [Steven's Julia intro](https://web.mit.edu/18.06/www/Fall17/1806/julia/Julia-intro.pdf))

\begin{align}V=\begin{bmatrix}1&\alpha _{1}&\alpha _{1}^{2}&\dots &\alpha _{1}^{n-1}\\1&\alpha _{2}&\alpha _{2}^{2}&\dots &\alpha _{2}^{n-1}\\1&\alpha _{3}&\alpha _{3}^{2}&\dots &\alpha _{3}^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&\alpha _{m}&\alpha _{m}^{2}&\dots &\alpha _{m}^{n-1}\end{bmatrix}\end{align}


```julia
using PyCall
np = pyimport("numpy")
```


```julia
np.vander(1:5, increasing=true)
```

[The source code for this function](https://github.com/numpy/numpy/blob/v1.16.1/numpy/lib/twodim_base.py#L475-L563) calls `np.multiply.accumulate` [which is implemented in C](https://github.com/numpy/numpy/blob/deea4983aedfa96905bbaee64e3d1de84144303f/numpy/core/src/umath/ufunc_object.c#L3678). However, this code doesn't actually perform the computation, it basically only checks types and stuff. The actual kernel is [implemented here](https://github.com/numpy/numpy/blob/deea4983aedfa96905bbaee64e3d1de84144303f/numpy/core/src/umath/loops.c.src#L1742), which is not even C code, but only a template that gets transformed to type-specific kernels. Still only a limited set of types `Float64`, `Float32`, and so forth are supported.

Here is a type-generic Julia implementation:


```julia
function vander(x::AbstractVector{T}) where T
    m = length(x)
    V = Matrix{T}(undef, m, m)
    for j = 1:m
        V[j,1] = one(x[j])
    end
    for i= 2:m
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
            end
        end
    return V
end
```


```julia
vander(1:5)
```

### Quick speed comparison


```julia
using BenchmarkTools
using Plots
ns = exp10.(range(1, 4, length=30));

tnp = Float64[]
tjl = Float64[]
for n in ns
    x = collect(1:n)
    push!(tnp, @belapsed np.vander($x) samples=3 evals=1)
    push!(tjl, @belapsed vander($x)    samples=3 evals=1)
end
plot(ns, tnp./tjl, m=:circle, xscale=:log10, yscale=:log10, ylims=[1, Inf],
     xlab="matrix size", ylab="NumPy time / Julia time", legend=:false)
```

Notably the clean and concise Julia implementation is **faster than numpy**'s C implementation for small matrices and **as fast** for large matrix sizes. Still it works for **arbitrary types**.


```julia
vander(Int32[4, 8, 16, 32])
```

This includes non-numerical types ... the only assumption is that the type induces a multiplicative group, i.e. has a `one` function to yield the identity element and an apropriate `*` defined.

A rather unusual one is `String`, which works since `one(String) == ""`.


```julia
vander(["this", "is", "a", "test"])
```

## How does this work?



 
- AST = Abstract Syntax Tree
- SSA = Static Single Assignment
- [LLVM](https://de.wikipedia.org/wiki/LLVM) = Low Level Virtual Machine

Julia compiles **ahead of time**, i.e. Julia waits until the first function call for a particular set of input types is issued and then compiles the full call stack all the way down to machine code specialising as much as possible along the way.


```julia
x = [1, 3, 5]  # Vector{Int64}

@time mysum(x)
@time mysum(x)
```

As we saw if thes change, Julia compiles a new specialization of the function!

To see what happens under the hood, we have a bunch of inspection macros:
* The AST after parsing (**`@macroexpand`**)
* The AST after lowering (**`@code_typed`**, **`@code_warntype`**)
* The AST after type inference and optimization (**`@code_lowered`**)
* The LLVM IR (**`@code_llvm`**)
* The assembly machine code (**`@code_native`**)


```julia
@code_typed mysum(1)
```


```julia
@code_lowered mysum(1)
```


```julia
@code_llvm mysum(1)
```

We can remove the comments (lines starting with `;` using `debuginfo=:none`).


```julia
@code_llvm debuginfo=:none mysum(1.2)
```


```julia
@code_native debuginfo=:none mysum(1.2)
```

Let's compare this to `Float64` input.


```julia
@code_native debuginfo=:none mysum(1)
```

## Types and code specialisation

As a rule of thumb: The Julia compiler only uses the types to specialise, not the values. Therefore using special types that match the possible assumptions most closely is crucial.

As an example we consider the determinant of a 2x2 matrix.


```julia
# Fast custom implementation ... note: No type annotation
det_custom(M) = @inbounds (M[1, 1] * M[2, 2] - M[1, 2] * M[2, 1])
```


```julia
M = randn(2, 2)
@btime det_custom(M);
```


```julia
using LinearAlgebra
@btime det(M);
```

Now that is a drastic difference ... but we actually did not tell the compiler anything about the size of `M`:


```julia
typeof(M)
```

Now let's help along ... by using a static array type:


```julia
using StaticArrays

S = @SMatrix randn(2, 2)
```


```julia
@btime det(S);
```

Got it! Now ... what happens under the hood:


```julia
@code_typed debuginfo=:none det(S)
```

in other words
```
%2  = S[1, 1]
%4  = S[2, 2]
%5  = S[1, 1] * S[2, 2]
%7  = S[2, 1]
%9  = S[1, 2]
%10 = S[2, 1] * S[1, 2]
%11 = S[1, 1] * S[2, 2] - S[2, 1] * S[1, 2]
```
exactly what we had hand-optimised before!

Now one should add here that part of the magic has not been the compiler being smart, but much rather the `StaticArrays` package, which has a number of carfully optimised implementations for standard operations:

```
============================================
    Benchmarks for 3×3 Float64 matrices
============================================
Matrix multiplication               -> 5.9x speedup
Matrix multiplication (mutating)    -> 1.8x speedup
Matrix addition                     -> 33.1x speedup
Matrix addition (mutating)          -> 2.5x speedup
Matrix determinant                  -> 112.9x speedup
Matrix inverse                      -> 67.8x speedup
Matrix symmetric eigendecomposition -> 25.0x speedup
Matrix Cholesky decomposition       -> 8.8x speedup
Matrix LU decomposition             -> 6.1x speedup
Matrix QR decomposition             -> 65.0x speedup
```

##### More details
- https://github.com/JuliaCI/BenchmarkTools.jl
- https://github.com/JuliaArrays/StaticArrays.jl/
- https://docs.julialang.org/en/v1/devdocs/compiler/

## Takeaways

- Julia can be fast
- Code is compiled down to the metal
- Stages can be inspected using macros like `@code_warntype`
- Typing is crucial to exploit specialised structure
- Type annotations are irrelevant for performance

0edaad4770c1e8932103abd6d77341079458b064
ipynb
Jupyter Notebook
05_Specialisation_Speed.ipynb
tuanvo-git/2022-rwth-julia-workshop
0b10ede0d49b4121adb9f6932a10623f7eef7015
MIT
05_Specialisation_Speed.ipynb
tuanvo-git/2022-rwth-julia-workshop
0b10ede0d49b4121adb9f6932a10623f7eef7015
MIT
05_Specialisation_Speed.ipynb
tuanvo-git/2022-rwth-julia-workshop
0b10ede0d49b4121adb9f6932a10623f7eef7015
MIT
Qwen/Qwen-72B
1. YES
2. YES