What to do when you are stuck with Julia?

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2024/01/12/rgg.html

Introduction

Today I decided to write about code refactoring in Julia.
This is a topic that is, in my experience, a quite big advantage of this language.

A common situation you are faced with when writing your code is as follows.
You need some functionality in your program and it is available in some library.
However, what the library provides does not meet your expectations.
Since in Julia most packages are written in Julia under MIT license, it is easy to solve this issue.
You just take the source code and modify it.

Today I want to show you a practical example of such a situation I have had this week
when working with the Graphs.jl package.

The post was written using Julia 1.10.0, BenchmarkTools.jl 1.4.0, and Graphs.jl 1.9.0.

The problem

In my work I needed to generate random geometric graphs.
This is a simple random graph model that works as follows
(here I describe a general idea, for details please check the
Wikipedia entry on random geometric graphs).
To generate a graph on N vertices you first drop N random points
in some metric space. Next you connect two points with an edge
if their distance is less than some pre-specified distance.

The Graphs.jl library provides the euclidean_graph function
that generates such graphs. Here is a summary of its docstring:

euclidean_graph(N, d; rng=nothing, seed=nothing, L=1., p=2., cutoff=-1., bc=:open)

Generate N uniformly distributed points in the box [0,L]^{d}
and return a Euclidean graph, a map containing the distance on each
edge and a matrix with the points' positions.

An edge between vertices x[i] and x[j] is inserted if norm(x[i]-x[j], p) < cutoff.
In case of negative cutoff instead every edge is inserted.
Set bc=:periodic to impose periodic boundary conditions in the box [0,L]^d.

So what is the problem with this function? Unfortunately it is slow.
Let us, for example check how long it takes to compute an average degree
of a node in such a graph with n nodes and cutoff=sqrt(10/n), when setting
bc=:periodic (periodic boundary, i.e. distance is measured on a torus) for two-dimensional space.

julia> using Graphs

julia> for n in 1_000:1_000:10_000
           println(@time ne(euclidean_graph(n, 2; cutoff=sqrt(10/n), bc=:periodic, seed=1)[1])/n)
       end
  0.091657 seconds (2.50 M allocations: 193.170 MiB, 14.06% gc time)
15.604
  0.300285 seconds (10.00 M allocations: 765.801 MiB, 12.59% gc time)
15.661
  0.686230 seconds (22.50 M allocations: 1.686 GiB, 12.47% gc time)
15.744
  1.175881 seconds (40.00 M allocations: 2.990 GiB, 10.97% gc time)
15.7065
  1.800561 seconds (62.50 M allocations: 4.666 GiB, 10.76% gc time)
15.6568
  2.697535 seconds (90.00 M allocations: 6.716 GiB, 14.76% gc time)
15.641333333333334
  3.690599 seconds (122.50 M allocations: 9.138 GiB, 13.28% gc time)
15.743571428571428
  4.745701 seconds (160.00 M allocations: 11.932 GiB, 13.53% gc time)
15.714
  5.962431 seconds (202.51 M allocations: 15.099 GiB, 12.70% gc time)
15.722222222222221
  7.195257 seconds (250.01 M allocations: 18.638 GiB, 11.42% gc time)
15.7086

We can see that the euclidean_graph function scales badly with n.
Note that by choosing cutoff=sqrt(10/n) we have a roughly constant
average degree, so the number of edges generates scales linearly, but
the generation time seems to grow much faster.

The investigation

To find out the source of the problem we can investigate the source
of euclidean_graph, which consists of two methods:

function euclidean_graph(
    N::Int,
    d::Int;
    L=1.0,
    rng::Union{Nothing,AbstractRNG}=nothing,
    seed::Union{Nothing,Integer}=nothing,
    kws...,
)
    rng = rng_from_rng_or_seed(rng, seed)
    points = rmul!(rand(rng, d, N), L)
    return (euclidean_graph(points; L=L, kws...)..., points)
end

function euclidean_graph(points::Matrix; L=1.0, p=2.0, cutoff=-1.0, bc=:open)
    d, N = size(points)
    weights = Dict{SimpleEdge{Int},Float64}()
    cutoff < 0.0 && (cutoff = typemax(Float64))
    if bc == :periodic
        maximum(points) > L && throw(
            DomainError(maximum(points), "Some points are outside the box of size $L")
        )
    end
    for i in 1:N
        for j in (i + 1):N
            if bc == :open
                Δ = points[:, i] - points[:, j]
            elseif bc == :periodic
                Δ = abs.(points[:, i] - points[:, j])
                Δ = min.(L .- Δ, Δ)
            else
                throw(ArgumentError("$bc is not a valid boundary condition"))
            end
            dist = norm(Δ, p)
            if dist < cutoff
                e = SimpleEdge(i, j)
                weights[e] = dist
            end
        end
    end
    g = Graphs.SimpleGraphs._SimpleGraphFromIterator(keys(weights), Int)
    if nv(g) < N
        add_vertices!(g, N - nv(g))
    end
    return g, weights
end

The beauty of Julia is that this source is written in Julia and is pretty short.
It immediately allows us to pinpoint the source of our problems. The core of we work
is done in a double-loop iterating over i and j indices. So the complexity of this
algorithm is quadratic in number of vertices.

Fixing the problem

The second beauty of Julia is that we can easily fix this. The idea can be found in
the Wikipedia entry on random geometric graphs in the algorithms section
here.

A simple way to improve the performance of the algorithm is to notice that if
you know L and cutoff you can partition the space into equal sized cells
having side length floor(Int, L / cutoff). Now you see that if you have a vertex
in some cell then it can be connected only to nodes in the same cell or cells directly
adjacent to it (the cells that are more far away contain the points that must be farther
than cutoff from our point). This means that we will have a much lower number of points
to consider. Below I show a code that is a modification of the original source that
adds this feature. The key added function is to_buckets which computes the bucket
identifier for each vertex and creates a dictionary that is a mapping from bucked
identifier to a vector of node numbers that fall into it:

using LinearAlgebra
using Random

function euclidean_graph2(
    N::Int,
    d::Int;
    L=1.0,
    rng::Union{Nothing,AbstractRNG}=nothing,
    seed::Union{Nothing,Integer}=nothing,
    kws...,
)
    rng = Graphs.rng_from_rng_or_seed(rng, seed)
    points = rmul!(rand(rng, d, N), L)
    return (euclidean_graph2(points; L=L, kws...)..., points)
end

function to_buckets(points::Matrix, L, cutoff)
    d, N = size(points)
    dimlen = max(floor(Int, L / max(cutoff, eps())), 1)
    buckets = Dict{Vector{Int}, Vector{Int}}()
    for (i, point) in enumerate(eachcol(points))
        bucket = floor.(Int, point .* dimlen ./ L)
        push!(get!(() -> Int[], buckets, bucket), i)
    end
    return buckets, dimlen
end

function euclidean_graph2(points::Matrix; L=1.0, p=2.0, cutoff=-1.0, bc=:open)
    d, N = size(points)
    weights = Dict{Graphs.SimpleEdge{Int},Float64}()
    cutoff < 0.0 && (cutoff = typemax(Float64))
    if bc == :periodic
        maximum(points) > L && throw(
            DomainError(maximum(points), "Some points are outside the box of size $L")
        )
    end
    buckets, dimlen = to_buckets(points, L, cutoff)
    deltas = collect(Iterators.product((-1:1 for _ in 1:size(points, 1))...))
    void = Int[]
    for (k1, v1) in pairs(buckets)
        for i in v1
            for d in deltas
                k2 = bc == :periodic ? mod.(k1 .+ d, dimlen) : k2 = k1 .+ d
                v2 = get(buckets, k2, void)
                for j in v2
                    i < j || continue
                    if bc == :open
                        Δ = points[:, i] - points[:, j]
                    elseif bc == :periodic
                        Δ = abs.(points[:, i] - points[:, j])
                        Δ = min.(L .- Δ, Δ)
                    else
                        throw(ArgumentError("$bc is not a valid boundary condition"))
                    end
                    dist = norm(Δ, p)
                    if dist < cutoff
                        e = Graphs.SimpleEdge(i, j)
                        weights[e] = dist
                    end
                end
            end
        end
    end

    g = Graphs.SimpleGraphs._SimpleGraphFromIterator(keys(weights), Int)
    if nv(g) < N
        add_vertices!(g, N - nv(g))
    end
    return g, weights
end

Note that it took less than 30 lines of code to add the requested feature to the code.

Performance of the improved code

Let us test our new code:

julia> for n in 1_000:1_000:10_000
           println(@time ne(euclidean_graph2(n, 2; cutoff=sqrt(10/n), bc=:periodic, seed=1)[1])/n)
       end
  0.017221 seconds (274.10 k allocations: 21.751 MiB, 22.15% gc time)
15.604
  0.022855 seconds (558.52 k allocations: 42.289 MiB, 10.43% gc time)
15.661
  0.032693 seconds (852.46 k allocations: 69.684 MiB, 8.52% gc time)
15.744
  0.043141 seconds (1.10 M allocations: 87.196 MiB, 14.73% gc time)
15.7065
  0.071273 seconds (1.41 M allocations: 109.725 MiB, 7.67% gc time)
15.6568
  0.068194 seconds (1.70 M allocations: 130.828 MiB, 12.54% gc time)
15.641333333333334
  0.071277 seconds (1.98 M allocations: 150.712 MiB, 11.85% gc time)
15.743571428571428
  0.081463 seconds (2.24 M allocations: 169.153 MiB, 10.67% gc time)
15.714
  0.099957 seconds (2.48 M allocations: 186.492 MiB, 8.08% gc time)
15.722222222222221
  0.148573 seconds (2.84 M allocations: 213.214 MiB, 18.37% gc time)
15.7086

We seem to get what we wanted. Our computation time looks to scale quite well.
Also the obtained average degree numbers are identical to the original ones.

Let us compare the performance on an even larger graph:

julia> n = 100_000;

julia> @time ne(euclidean_graph(n, 2; cutoff=sqrt(10/n), bc=:periodic, seed=1)[1])/n
908.976252 seconds (25.00 G allocations: 1.819 TiB, 11.64% gc time, 0.00% compilation time)
15.70797

julia> julia> @time ne(euclidean_graph2(n, 2; cutoff=sqrt(10/n), bc=:periodic, seed=1)[1])/n
  1.640495 seconds (27.53 M allocations: 2.121 GiB, 19.83% gc time)
15.70797

Indeed we see that the timing of the original implementation becomes prohibitive for
larger graphs.

Making sure things are correct

Before we finish there is one important task we need to make. We should check if
indeed the euclidean_graph2 function produces the same results as euclidean_graph.
This is easy to test with the following randomized test:

julia> using Test

julia> Random.seed!(1234);

julia> @time for i in 1:1000
           N = rand(10:500)
           d = rand(1:5)
           L = rand()
           p = 3*rand()
           cutoff = rand() * L / 4
           bc = rand([:open, :periodic])
           seed = rand(UInt32)
           @test euclidean_graph(N, d; L, p, cutoff, bc, seed) ==
                 euclidean_graph2(N, d; L, p, cutoff, bc, seed)
       end
 16.955773 seconds (275.09 M allocations: 20.342 GiB, 12.27% gc time)

We have tested 1000 random setups of the experiments. In each of them both functions
returned the same results.

Conclusions

In my post I have shown you an example that one can easily tweak some package code to your needs.
In this case this was performance, but it can be equally well functionality.

I did not comment too much on the code itself, as it was a bit longer than usual, but let me
discuss as a closing remark one performance aspect of the code. In my to_buckets function I used
the get! function to populate the dictionary with a mutable default value (Int[] in that case).
You might wonder why I preferred to use an anonymous function instead of passing a default
as a third argument. The reason is number of allocations. Check this code:

julia> using BenchmarkTools

julia> function f1()
           d = Dict(1 => Int[])
           for i in 1:10^6
               get!(d, 1, Int[])
           end
           return d
       end
f1 (generic function with 1 method)

julia> function f2()
           d = Dict(1 => Int[])
           for i in 1:10^6
               get!(() -> Int[], d, 1)
           end
           return d
       end
f2 (generic function with 1 method)

julia> @benchmark f1()
BenchmarkTools.Trial: 195 samples with 1 evaluation.
 Range (min … max):  19.961 ms … 45.328 ms  ┊ GC (min … max): 10.26% … 13.19%
 Time  (median):     23.962 ms              ┊ GC (median):    10.45%
 Time  (mean ± σ):   25.660 ms ±  5.050 ms  ┊ GC (mean ± σ):  12.27% ±  3.76%

    ▃▂▂█▃▃ ▃▂▃
  ▆▄██████▇███▇▆▄▄▇▄▄▁▃▆▄▃▅▄▅▄▄▄▃▄▃▁▃▃▁▁▁▃▃▃▃▃▃▁▃▁▁▃▁▁▁▁▁▃▁▁▃ ▃
  20 ms           Histogram: frequency by time        44.3 ms <

 Memory estimate: 61.04 MiB, allocs estimate: 1000005.

julia> @benchmark f2()
BenchmarkTools.Trial: 902 samples with 1 evaluation.
 Range (min … max):  4.564 ms … 11.178 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     5.149 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.526 ms ±  1.396 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▅▄▄▃▃▄▅▆▄▁▁                                           ▃
  ██████████████▆▆▆▆▆▅▇▆▄▆▆▅▅▁▄▅▁▅▁▇▅▆▄▄▇▁▅▅▅▄▁▄▄▆▅▁▅▁▁▅▆██▅ █
  4.56 ms      Histogram: log(frequency) by time       10 ms <

 Memory estimate: 592 bytes, allocs estimate: 5.

As you can see f2 is much faster than f1 and does much less allocations. The issue
is that f1 allocates a fresh Int[] object in every iteration of the loop, while
f2 would allocate it only if get! does not hit a key that already exists in d
(and in our experiment I always queried for 1 which was present in d).

Code cleanup for the New Year

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2024/01/05/mastermind.html

Introduction

I have done a small cleanup of my folders with code for the New Year.
In the process I have discovered some old Julia code that I have written and forgotten about.
One of such scripts is a simple Mastermind solver.
So I decided to share it in this post.

The code was tested under Julia 1.10.0, Combinatorics.jl 1.0.2, FreqTables.jl 0.4.6, and StatsBase.jl 0.34.2.

The Mastermind game

I assume you know how Mastermind is played.
If not – I recommend you read the Wikipedia description before continuing.

Let me briefly summarize the version of the game that I used as a reference
(this is a version that was popular in Poland when I was young):

  • The game is played using 8 colors of pegs.
  • The code that we are required to break (a solution) consists of 5 pegs of which each has to have a different color.
  • You are allowed to guess the solution by showing 5 pegs (that also have to be of different colors)
    and you get two numbers: black representing number of pegs matching color and position between your guess and a solution,
    and white dots representing pegs that have a correct color, but wrong position.
  • The game is finished when your guess matches the solution (i.e. you get 5 blacks and 0 whites).
    The outcome of the game is the number of guesses you have made.

The problem

In my code I wanted to compare two strategies. Both of them use a concept of available options,
so let us define it first.

First note that we have 5 pegs and each of them has one of 8 unique colors that we have:

julia> binomial(8, 5) * factorial(5)
6720

initially available options as solution when we start a game. After each move we can check which of options
match the responses to your guesses. This matching reduces the number of available options.

The first strategy we will consider is to randomly pick a guess that is still an available option.
This strategy is in practice possible to follow (up to ensuring randomness of choice which can be hard :))
even without the help of the computer. At any moment during the game you just need to show a guess that
is consistent with the answers you already got.

The second strategy, called minimax, works as follows. Before making a move we consider every initially
available option as a possible candidate. For each candidate we check what black-white combination
we would get against all available options. From this we can get a distribution of number of black-white
combinations that a given guess would produce. We sort this distribution in a descending order and pick
the guess that produces the lexicographically smallest value.

Let me explain this by example. Assume we consider guesses A, B, and C and have still 20 available options.

For the guess A we count that we could have 10 answers of the form 2 black & 2 white, 5 answers 1 black & 2 white,
and 5 answers 2 black & 1 white. Thus for option A the encoding is [10, 5, 5].

For the guess B we count that we could have 3 answers of the form 2 black & 3 white, 5 answers 1 black & 2 white,
and 12 answers 2 black & 1 white. Thus for option B the encoding is [12, 5, 3].

For the guess C we count that we could have 10 answers of the form 1 black & 1 white, 8 answers 2 black & 2 white,
and 2 answers 1 black & 1 white. Thus for option C the encoding is [10, 8, 2].

Therefore since lexicographically [10, 5, 5] < [10, 8, 2] < [12, 5, 3] the option A is best (and we would pick it)
and the option B is worst.

The rationale behind this strategy is that we want to minimize the maximum possible number of available options we will be
left with in the next round.

Clearly minimax strategy is much more complex. For sure it is hard to be used without computer assistance. The question is
if it is better than random strategy and, if so, by how much.

The code

Let us start with the code that implements the game.

The first part is environment setup:

using Combinatorics
using FreqTables
using Random
using StatsBase

struct Position
    pos::NTuple{5, Int8}
    pegs::NTuple{8, Bool}
end

const ALL_POSITIONS = let
    combs = combinations(Int8.(1:8), 5)
    perms = reduce(vcat, (collect∘permutations).(combs))
    [Position(Tuple(pos), ntuple(i -> i in pos, 8)) for pos in perms]
end

In this part of code we loaded all required packages. Next we define Position structure.
It encodes the colors (as values from 1 to 8) passed for each peg. The pos field
directly encodes the position. The pegs field is a helper code containing information
if a given color was used in pos (to make it simpler to later use this data structure).
Finally we create ALL_POSITIONS vector which contains all 6720 available options at the start of the game.

Next define a function counting the number of black and white matches between the guess and the solution:

function match(p1::Position, p2::Position)
    pos1, pegs1 = p1.pos, p1.pegs
    pos2, pegs2 = p2.pos, p2.pegs
    black = sum(x -> x[1] == x[2], zip(p1.pos, p2.pos))
    matches = sum(x -> x[1] & x[2], zip(p1.pegs, p2.pegs))
    return black, matches - black
end

Here we see why both pos and pegs values are useful. We use pos to compute the number of black matches.
Then we use pegs to compute the number of all matches (disregarding the location). Therefore matches - black
gives ut the number of white matches.

We are now ready to implement our decision rule functions:

random_rule(options::AbstractVector{Position}) = rand(options)

function minimax_rule(options::AbstractVector{Position})
    length(options) == 1 && return only(options)
    best_guess = options[1]
    best_guess_lengths = [length(options)]
    for guess in shuffle(options)
        states = Dict{Tuple{Int,Int},Int}()
        for opt in options
            resp = match(guess, opt)
            states[resp] = get(states, resp, 0) + 1
        end
        current_guess_lengths = sort!(collect(values(states)), rev=true)
        if current_guess_lengths < best_guess_lengths
            best_guess = guess
            best_guess_lengths = current_guess_lengths
        end
    end
    return best_guess
end

The random_rule is simple. We just randomly pick one of the available options.

The minimax_rule code is more complex. We traverse guess taken from options in random order.
For each guess in the states dictionary we keep the count of black-white answers generated by
this guess. Then current_guess_lengths keeps the vector of combination counts sorted in a descending order.
After checking all guess values we return the best_guess value which holds the guess with the smallest value of
current_guess_lengths.

The final part of the code is a rule testing function:

function run_test(rule::Function, solution::Position)
    options = ALL_POSITIONS
    moves = 0
    while true
        moves += 1
        guess::Position = rule(options)
        guess == solution && return moves
        resp = match(guess, solution)
        options = filter(opt -> match(guess, opt) == resp, options)
        @assert solution in options
    end
end

In the run_test we get the rule function and the solution to the game. Next iteratively we ask the rule function what next move it picks. If it is the solution
then we return the number of moves taken. Otherwise we filter the available options, stored in the options variable to leave only those that match the responses
received by matching the guess against the solution.

Now the big question is if the extra complexity of minimax_rule is worth the effort?

The test

Now we are ready to test the random_rule against the minimax_rule. One important observation that we can make is that it is enough to run the run_test function
with only one starting solution, which I pick to be ALL_OPTIONS[1]. The reason is that the problem is fully symmetric (i.e. all possible solutions are indistinguishable with respect to recoloring)
and both random_rule and minimax_rule are fully randomized. For random_rule it is clear. For minimax_rule the key line of code that ensures it is for guess in shuffle(options).
In this way we ensure that we do not favor any options. E.g. if we did for guess in options then clearly we would favor options[1] option (as the game would always finish in the first move).

Now let us run the tests:

julia> Random.seed!(1234);

julia> @time random_res = [run_test(random_rule, ALL_POSITIONS[1]) for i in 1:1_000_000];
138.395938 seconds (7.13 M allocations: 97.072 GiB, 1.59% gc time, 0.11% compilation time)

julia> @time minimax_res = [run_test(minimax_rule, ALL_POSITIONS[1]) for i in 1:1_000];
1162.463279 seconds (65.08 M allocations: 20.119 GiB, 0.03% gc time, 0.04% compilation time)

As you can see I run one thousand times more random_rule than minimax_rule as the former is much faster (and still ten times more time is spent in gathering the results for minimax_rule).

Now let us collect some statistics about the both approaches to the solution:

julia> proptable(random_res)
9-element Named Vector{Float64}
Dim1  │
──────┼─────────
1     │ 0.000156
2     │ 0.002364
3     │  0.02251
4     │ 0.139851
5     │  0.39269
6     │ 0.350372
7     │  0.08736
8     │ 0.004646
9     │   5.1e-5

julia> proptable(minimax_res)
6-element Named Vector{Float64}
Dim1  │
──────┼──────
3     │ 0.015
4     │ 0.138
5     │ 0.476
6     │ 0.336
7     │ 0.034
8     │ 0.001

Note that the probability of solving the puzzle in one move is (we hit the solution by chance):

julia> 1 / (binomial(8, 5) * factorial(5))
0.00014880952380952382

It is pretty accurately approximated for random_res. For minimax_res we did not get even one such observation because our sample was too small.
In general, what we can see is that the probabilities of numbers of moves are similar (keeping in mind the fact that minimax_res has only 1000 observations),
except that random rule less frequently finishes in 5 moves, and more frequently in 7 moves. However, the difference is not very large.
This is also reflected by the fact that the mean time required to finish the game is only slightly lower for minimax strategy:

julia> mean(random_res)
5.346647

julia> mean(minimax_res)
5.239

Conclusions

Interestingly, for the variant of Mastermind we considered a random strategy is not much worse than a much more expensive minimax approach.
A limitation of our minimax approach is that it is doing only one-step-ahead optimization. As an extension to our exercise you could consider
checking how much better would be a full-depth minimax algorithm.

Happy hacking!

Personal Julia 2023 retrospective

By: julia on Abel Soares Siqueira

Re-posted from: https://abelsiqueira.com/blog/2024-01-04-personal-julia-2023-retrospective/

2023 was a very interest year for me, regarding the use and promotion of the Julia Programming Language. Here are the three main reasons why.
First, I gave no less than 4 talks (not counting anything in-company).
On 19/January at the JuliaLang Eindhoven Meetup. I talked about the difficulties of maintaining a large ecosystems of packages. It was the first time I went to Eindhoven, and I met some great people there.