Generic programming¶

From Wikipedia:

Generic programming is a style of computer programming in which algorithms are written in terms of types to-be-specified-later that are then instantiated when needed for specific types provided as parameters.

Example: Vandermonde matrix (revisited)¶

Vandermonde matrix: \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}

Keep function argument types generic if possible¶

In [ ]:

function vander_naive(x::Vector)
    m = length(x)
    V = Matrix{Float64}(undef, m, m)
    for j = 1:m
        V[j,1] = 1.0
    end
    for i= 2:m
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
            end
        end
    return V
end

In [ ]:

x = rand(3)
vander_naive(x)

In [ ]:

vander_naive(1:3)

The annotation ::Vector in the the function signature is unnecessarily specific.

In [ ]:

function vander_naive(x::AbstractVector)
    m = length(x)
    V = Matrix{Float64}(undef, m, m)
    for j = 1:m
        V[j,1] = 1.0
    end
    for i= 2:m
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
            end
        end
    return V
end

In [ ]:

vander_naive(1:3)

Avoid explicit typing if possible¶

In [ ]:

vander_naive([1,2,3])

Why is the result a matrix of floating point numbers....?

Even worse:

In [ ]:

vander_naive(rand(ComplexF64, 3))

We can easily cover those cases as well by only slightly modifying our code.

In [ ]:

function vander_almost_generic(x::AbstractVector)
    T = eltype(x)
    m = length(x)
    V = Matrix{T}(undef, m, m)
    for j = 1:m
        V[j,1] = 1.0
    end
    for i= 2:m
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
            end
        end
    return V
end

In [ ]:

vander_almost_generic([1,2,3])

In [ ]:

vander_almost_generic(rand(ComplexF64, 3))

In [ ]:

vander_almost_generic(["Stadt", "Land", "Fluss"])

In [ ]:

function vander_generic(x::AbstractVector{T}) where T # this is the same as just x::AbstractVector
    m = length(x)
    V = Matrix{T}(undef, m, m)
    for j = 1:m
        V[j,1] = one(T)
    end
    for i= 2:m
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
            end
        end
    return V
end

In [ ]:

vander_generic(["Stadt", "Land", "Fluss"])

One more level of generality, just because we can. :)

In [ ]:

vander_generic([3, "Stadt", 4 + 5im])

In [ ]:

function vander_supergeneric(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

In [ ]:

vander_supergeneric([3, "Stadt", 4 + 5im])

And all of this comes at no performance penality¶

In [ ]:

using BenchmarkTools

In [ ]:

x = rand(Float64, 100);
@btime vander_naive($x);
@btime vander_supergeneric($x);

Actually, for this specific example our generic code is faster in a few cases inasmuch as type conversions are unnecessary.

In [ ]:

x = rand(Int, 100);
@btime vander_naive($x);
@btime vander_supergeneric($x);

In [ ]:

x = rand(Bool, 100);
@btime vander_naive($x);
@btime vander_supergeneric($x);

On the other hand, sometimes it is worth converting to a different type to dispatch to a faster method or to utilize magic like compiler optimizations.

In [ ]:

x = rand(Float32, 100);
@btime vander_naive($x);
@btime vander_supergeneric($x);

Arbitrary precision computations¶

Let's say you have implemented the following crazily complex physics code as part of your thesis project which takes in a number and spits out the answer to life, the universe, and everything.

In [ ]:

function answer_to_life_universe_and_everything(x::Integer)
    m = sin((2*x)^100)
    c = 42
    E = m*c^2
    a = sqrt(abs(E))
    b = atan(m)
    c = a^2 + b^2
    answer = sqrt(1764)/(1+exp(-c))
end

In [ ]:

answer_to_life_universe_and_everything(2)

Alright, apparently the answer is 21!

The author: "I checked the code multiple times, it is correct. So let's publish."

Without changing a line of code we can check the correctness (in the numerical sense) of our result using arbitrary precision arithmetics.

In [ ]:

big(2)

In [ ]:

typeof(big(2))

In [ ]:

answer_to_life_universe_and_everything(big(2))

Let's get a bit more fancy!¶

In [ ]:

using Interact

In [ ]:

@manipulate for n in 1:20
    [i*j for i in 1:n, j in 1:n]
end

In [ ]:

function insert_block(A::AbstractMatrix, i, j, what=7)
    B = copy(A)
    B[i:i+2, j:j+2] .= what
    B
end

In [ ]:

A = fill(0, 9, 9)
insert_block(A, 3, 5)  # this returns the new matrix

In [ ]:

A = fill(0, 10, 10)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j)
end

Let's add some color!¶

Our function insert_block is generic. Since the first argument A isa AbstractArray, we can index into it and set new values. Pretty much every value type is fine!

In [ ]:

using Colors

In [ ]:

@manipulate for n in 1:80
    distinguishable_colors(n)
end

In [ ]:

colors = distinguishable_colors(10)

In [ ]:

colors[1]

In [ ]:

A = fill(colors[1], 10, 10)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j, colors[4])
end

Generic Programming + Multiple Dispatch + JIT¶

The possibility to write generic algorithms that compile to fast machine code in combination with multiple dispatch leads to an (unreasonable) amount of code reuse. This sharing of code comes in two forms:

Sharing types among a wide variety of packages implementing different algorithms;
Sharing generic algorithms that work for different package-defined types implementing common abstractions.

Sharing types: DataStructures.jl, OrderedCollections.jl, StaticArrays.jl, Colors.jl, Measurements.jl ...

As of the time of this writing, 804 packages depend on the data types provided in DataStructures.jl.

851 packages reuse type implementations in OrderedCollections.jl.

That's about every third package!

Sharing generic algorithms: StatsBase.jl, SortingAlgorithms.jl, GenericLinearAlgebra.jl, ...

Emergent features¶

Measurement type from Measurements.jl and differential equation solver from OrdinaryDiffEq.jl (i.e. DifferentialEquations.jl)

In [ ]:

using OrdinaryDiffEq, Measurements, PyPlot

#Half-life of Carbon-14 is 5730 years.
c = 5.730 ± 2

#Setup
u0 = 1.0 ± 0.1
tspan = (0.0 ± 0.0, 1.0 ± 0.0)

#Define the problem
radioactivedecay(u,p,t) = -c*u

#Pass to solver
prob = ODEProblem(radioactivedecay,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

In [ ]:

# analytic solution
u = u0 .* exp.(-c .* sol.t);

# plot solution
ts = getfield.(sol.t, :val)
solvals = getfield.(sol, :val)
solerrs = getfield.(sol, :err);

errorbar(ts, solvals, yerr=solerrs)
plot(ts, getfield.(u, :val), color="red", lw=2)
ylabel("u(t)")
xlabel("t");

Note that, in some sense, Julia implemented that feature by itself.

The authors of Measurements.jl and DifferentialEquations.jl never had any collabration on this.

It just works.

Core messages of this Notebook¶

It is simple to write type-generic code in Julia and you should do it.
Generally, generic code is just as fast as specific code.
Generic Programming + Multiple Dispatch + JIT = lots of code sharing and emergent features