By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2023/11/10/pe215.html

Introduction

Recently I have written several posts on tips & tricks regarding the usage of Base Julia functionalities.
Today I thought of presenting a solution to a bigger problem that would feature some of the basic data structures available in Julia.

For this I have chosen to solve the Project Euler problem 215.

The post was written under Julia 1.9.2, Graphs.jl 1.9.0, and Chain.jl 0.5.0.

The problem

The Project Euler problem 215 is stated as follows:

Consider the problem of building a wall out of 2×1 and 3×1 bricks (horizontal x vertical dimensions)
such that, for extra strength, the gaps between horizontally-adjacent bricks never line up in consecutive layers,
i.e. never form a “running crack”.
There are eight ways of forming a crack-free 9×3 wall, written W(9, 3) = 8. Calculate W(32, 10).

(I encourage you to visit the source web page for a visualization.)

The solution approach

The basic approach to the solution of the problem can be done in the following steps:

1. Compute the list of ways a single layer of the wall can be built.
2. Compute which layers can be adjacent in the wall.
3. Compute the number of possibilities of building the full wall.

Let us go through these steps on the 9×3 case for which we know the solution.

Compute the list of ways a single layer can be built

For the first step we define the following function:

function list_valid_layers(n::Integer)
    function build(part, len)
        if n-1 <= len + 2 <= n
            push!(valid, BitSet(part))
        end
        for i in 2:3
            if len + i < n
                push!(part, len + i)
                build(part, len + i)
                pop!(part)
            end
        end
    end

    valid = BitSet[]
    build(Int[], 0)
    return valid
end

It takes the layer length n as an input and returns a vector of possible compositions of a layer.
Note that in order to identify the composition of a layer it is enough to keep track of where
the “cracks” in the wall are. For performance we keep this information in BitSet.

The list of compositions is created recursively in the build inner functions. We first check if the
wall would be finished by adding a brick of width 2 or 3. If yes, we store a BitSet pattern of cracks
in the valid vector. If not we try to add another brick of width 2 or 3 to the wall and call the
build function recursively. Note the use of push! and pop! functions in this part of code. It allows
us to reuse a single part vector to keep track of the state of the computations.

Let us check the list of valid layers of width 9:

julia> v9 = list_valid_layers(9)
5-element Vector{BitSet}:
 BitSet([2, 4, 6])
 BitSet([2, 4, 7])
 BitSet([2, 5, 7])
 BitSet([3, 5, 7])
 BitSet([3, 6])

We see that we have 5 such ways. In one of them there are 3 bricks of with 3.
In the remaining we have a single brick of with 3 and four bricks of width 2 (giving 4 ways to build a layer).

Compute which layers can be adjacent in the wall

Two layers can be adjacent in our wall if the intersection of their cracks is empty. Let us keep track of this information
in a graph:

julia> using Graphs

julia> function build_graph(valid::Vector)
           g = Graph(length(valid))
           for i in 1:length(valid), j in i+1:length(valid)
               isempty(valid[i] ∩ valid[j]) && add_edge!(g, i, j)
           end
           return g
       end
build_graph (generic function with 1 method)

julia>

julia> g9 = build_graph(v9)
{5, 3} undirected simple Int64 graph

julia> collect(edges(g9))
3-element Vector{Graphs.SimpleGraphs.SimpleEdge{Int64}}:
 Edge 1 => 4
 Edge 2 => 5
 Edge 3 => 5

We see that there are only three allowed pairs of adjacent layers:

• [2, 4, 6] and [3, 5, 7];
• [2, 4, 7] and [3, 6];
• [2, 5, 7] and [3, 6].

Note that when building a graph we used a simple graph. The reason is that the adjacency relation is symmetric
(if layer i can follow layer j then the opposite is also true).

Finally note that I iterate i and j assuming 1-based indexing of valid. But this assumption is met as I
require that it is a Vector.

Compute the number of possibilities of building the full wall

What is left is computing the final number of allowed layer compositions.
Here is a function that performs this task:

function count_layers(graph::SimpleGraph, depth::Integer)
    layer_count = fill(1, nv(graph), depth)
    for k in 2:depth, i in 1:nv(graph)
        layer_count[i, k] = sum(layer_count[j, k-1] for j in neighbors(graph, i))
    end
    return layer_count
end

Note that we keep the result in a layer_count matrix. Its layer_count[i, k] entry
tells us in how many ways we can construct a wall having k layers that ends with layer i
(where i is the location of layer in the vector returned by list_valid_layers).

To get the desired result we will need to count the sum of all possibilities from the final layer:

julia> count_layers(g9, 3)
5×3 Matrix{Int64}:
 1  1  1
 1  1  2
 1  1  2
 1  1  1
 1  2  2

julia> sum(count_layers(g9, 3)[:, end])
8

We can see that the result is 8, which is the same as in the problem statement.
So, hopefully we are ready to solve the problem of computing W(32, 10).

Computation of W(32, 10)

To compute W(32, 10) we need to invoke the functions we have defined in this post in sequence.
This kind of task is a natural application of Chain.jl @chain function. So let us test it:

using Chain

@chain 32 begin
    list_valid_layers
    build_graph
    count_layers(10)
    sum(_[:, end])
end

As usual in this blog, I am not presenting the answer, to encourage you to try solving the puzzle yourself.
However, let me remark that the obtained number is relatively large, so you should carefully check if we do not hit the Int64 overflow issue.
If you are not sure what the problem could be check out my recent post, where I explain it.

Conclusions

As a conclusion let us check how long it takes to solve the problem on my laptop:

julia> @elapsed @chain 32 begin
           list_valid_layers
           build_graph
           count_layers(10)
           sum(_[:, end])
       end
0.4039108

The timing of around 0.4 seconds is not that bad, bearing in mind that we used a brute force way to solve the problem.
However, to get this timing keep in mind the following data structures that we used
(not all of them were performance bottleneck in this problem, but I think it is worth to recap them):

• push! and pop! updates of a vector in a recursion (limiting the number of allocations we needed to make);
• use of BitSet to keep track of the “cracks” in a layer (taking advantage of the fact that our data was small);
• use of SimpleGraph to efficiently navigate the allowed links between walls;
• use of a matrix to incrementally count the number of allowed options.

In all of these steps we took advantage of the fact that for loops in Julia are fast. Note that I took care
to put the code inside functions not only to have reusable components, but also to make sure that
Julia will efficiently execute the code.

Enjoy!

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2023/11/03/vec.html

Introduction

Last week I have written a “Learning to zip stuff in Julia” post
that discussed the zip function. To my surprise, even though the topic was basic,
it has received a lot of positive feedback. Therefore I thought of writing about
another entry-level problem.

Often, when working with arrays in Julia we want to flatten them to a vector.
Today, I want to discuss two ways how you can do it and the differences between them.

The post was written under Julia 1.9.2.

The problem

Assume you have the following three dimensional array:

julia> a3d = [(i, j, k) for i in 1:2, j in 1:3, k in 1:4]
2×3×4 Array{Tuple{Int64, Int64, Int64}, 3}:
[:, :, 1] =
 (1, 1, 1)  (1, 2, 1)  (1, 3, 1)
 (2, 1, 1)  (2, 2, 1)  (2, 3, 1)

[:, :, 2] =
 (1, 1, 2)  (1, 2, 2)  (1, 3, 2)
 (2, 1, 2)  (2, 2, 2)  (2, 3, 2)

[:, :, 3] =
 (1, 1, 3)  (1, 2, 3)  (1, 3, 3)
 (2, 1, 3)  (2, 2, 3)  (2, 3, 3)

[:, :, 4] =
 (1, 1, 4)  (1, 2, 4)  (1, 3, 4)
 (2, 1, 4)  (2, 2, 4)  (2, 3, 4)

In many situations you might want to transform it into a one dimensional vector.
For example many functions explicitly require AbstractVector as their input.
The question is how can you do it.

There are two fundamental ways to perform this operation. The first one creates
a new independent vector from the source data, and the second one reuses the memory
of the source data. Let me discuss them in more detail.

Copied vector

If you want to copy the data in a3d into a new vector you can write:

julia> a3d[:]
24-element Vector{Tuple{Int64, Int64, Int64}}:
 (1, 1, 1)
 (2, 1, 1)
 (1, 2, 1)
 (2, 2, 1)
 (1, 3, 1)
 (2, 3, 1)
 (1, 1, 2)
 (2, 1, 2)
 (1, 2, 2)
 (2, 2, 2)
 (1, 3, 2)
 (2, 3, 2)
 (1, 1, 3)
 (2, 1, 3)
 (1, 2, 3)
 (2, 2, 3)
 (1, 3, 3)
 (2, 3, 3)
 (1, 1, 4)
 (2, 1, 4)
 (1, 2, 4)
 (2, 2, 4)
 (1, 3, 4)
 (2, 3, 4)

The syntax is short and easy to read. It takes advantage of the fact that every Julia
array should support linear indexing, as is explained in the Julia Manual section
on Linear Indexing.

The benefit of this approach is that the new object is freshly allocated, so modifying
it will not modify the source. The downside is that it requires memory allocation.

Aliased vector

If we want to avoid excessive memory allocation (which might be relevant for large objects)
we can create a vector from an array without copying the data. This can be achieved using
the vec function:

julia> v = vec(a3d)
24-element Vector{Tuple{Int64, Int64, Int64}}:
 (1, 1, 1)
 (2, 1, 1)
 (1, 2, 1)
 (2, 2, 1)
 (1, 3, 1)
 (2, 3, 1)
 (1, 1, 2)
 (2, 1, 2)
 (1, 2, 2)
 (2, 2, 2)
 (1, 3, 2)
 (2, 3, 2)
 (1, 1, 3)
 (2, 1, 3)
 (1, 2, 3)
 (2, 2, 3)
 (1, 3, 3)
 (2, 3, 3)
 (1, 1, 4)
 (2, 1, 4)
 (1, 2, 4)
 (2, 2, 4)
 (1, 3, 4)
 (2, 3, 4)

Now v is a vector that shares memory with a3d. We can check it using the pointer function:

julia> pointer(v)
Ptr{Tuple{Int64, Int64, Int64}} @0x000002606d121540

julia> pointer(a3d)
Ptr{Tuple{Int64, Int64, Int64}} @0x000002606d121540

or the Base.mightalias function:

julia> Base.mightalias(a3d, v)
true

The benefit is that vec(a3d) is in general much faster than a3d[:] since it makes less work.
The downside is that you need to remember that mutating one of the objects will change the other.
For example:

julia> a3d[1, 1, 1] = (100, 100, 100);

julia> a3d
2×3×4 Array{Tuple{Int64, Int64, Int64}, 3}:
[:, :, 1] =
 (100, 100, 100)  (1, 2, 1)  (1, 3, 1)
 (2, 1, 1)        (2, 2, 1)  (2, 3, 1)

[:, :, 2] =
 (1, 1, 2)  (1, 2, 2)  (1, 3, 2)
 (2, 1, 2)  (2, 2, 2)  (2, 3, 2)

[:, :, 3] =
 (1, 1, 3)  (1, 2, 3)  (1, 3, 3)
 (2, 1, 3)  (2, 2, 3)  (2, 3, 3)

[:, :, 4] =
 (1, 1, 4)  (1, 2, 4)  (1, 3, 4)
 (2, 1, 4)  (2, 2, 4)  (2, 3, 4)

julia> v
24-element Vector{Tuple{Int64, Int64, Int64}}:
 (100, 100, 100)
 (2, 1, 1)
 (1, 2, 1)
 (2, 2, 1)
 (1, 3, 1)
 (2, 3, 1)
 (1, 1, 2)
 (2, 1, 2)
 (1, 2, 2)
 (2, 2, 2)
 (1, 3, 2)
 (2, 3, 2)
 (1, 1, 3)
 (2, 1, 3)
 (1, 2, 3)
 (2, 2, 3)
 (1, 3, 3)
 (2, 3, 3)
 (1, 1, 4)
 (2, 1, 4)
 (1, 2, 4)
 (2, 2, 4)
 (1, 3, 4)
 (2, 3, 4)

Conclusions

In summary if in your analytical workflow you need a vector,
but you have a general array a as an input you have two options:

• write a[:], which creates a copy (safer, but slower and uses more memory);
• write vec(a), which creates an alias (unsafe, but faster and uses less memory).

In practice I find both options useful depending on the circumstances, so I think it is worth to be aware of them.

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2023/10/27/zip.html

Introduction

Today I want to discuss the basics of the design of the zip function in Julia.
This is a relatively introductory topic but, in my experience, wrong usage of this function can lead to hard to catch bugs.

The post was written under Julia 1.9.2.

What does zip do?

Let us start with the description of the zip function from its documentation:

Run multiple iterators at the same time, until any of them is exhausted.
The value type of the zip iterator is a tuple of values of its subiterators.

If it is not fully clear what it does let me give two examples:

julia> [v for v in zip(1:3, 11:13)]
3-element Vector{Tuple{Int64, Int64}}:
 (1, 11)
 (2, 12)
 (3, 13)

julia> [v for v in zip("abcdef", Iterators.cycle([1, 2]))]
6-element Vector{Tuple{Char, Int64}}:
 ('a', 1)
 ('b', 2)
 ('c', 1)
 ('d', 2)
 ('e', 1)
 ('f', 2)

The examples show the following nice features of the zip function:

• It can work with any iterator, not just e.g. vectors.
  In the second example we used a string and a special cyclic iterator from Iterators module.
  Importantly we can work with iterators that are not indexable.
• The zip function works until the shortest iterator of all passed to it is exhausted.
  It is often useful when we have some of the iterators that are infinite.
  Again, in our second example Iterators.cycle([1, 2])) is infinite and allowed us to clearly indicate even and odd characters in the string.

The last nice feature of zip is that it is lazy. To understand what it means check this:

julia> zip(1:3, 11:13)
zip(1:3, 11:13)

julia> typeof(zip(1:3, 11:13))
Base.Iterators.Zip{Tuple{UnitRange{Int64}, UnitRange{Int64}}}

We can see that the zip function by itself does not do any work. Instead, it returns an iterator that can be queried.
That is why in the first examples I used comprehensions to show the values produced by the iterator returned by zip.
Why is it useful? Being lazy avoids excessive memory allocation, which is often important, especially if you work with large data.

Why is zip risky?

So what are the potential problems with zip in practice? The issue is that we often want to use it in scenarios
when we would want to be sure that all iterators passed to zip have the same length. Unfortunately zip does not check it.

Here is a simple example of a problematic definition of a function comupting a dot product:

julia> dot_bad(x, y) = sum(a * b for (a, b) in zip(x, y))
dot_bad (generic function with 1 method)

julia> dot_bad(1:3, 4:6)
32

julia> dot_bad(1:3, 4:7)
32

julia> dot_bad(1:4, 4:6)
32

Typically with such a function we would want it to error if it gets input iterators of unequal length.

In practice, if it is required that the passed iterators must have equal length, they usually support the length function.
Therefore often fix of the definition that is similar to the following is enough:

julia> function dot_safer(x, y)
           @assert length(x) == length(y)
           return sum(a * b for (a, b) in zip(x, y))
       end
dot_safer (generic function with 1 method)

julia> dot_safer(1:3, 4:6)
32

julia> dot_safer(1:3, 4:7)
ERROR: AssertionError: length(x) == length(y)

julia> dot_safer(1:4, 4:6)
ERROR: AssertionError: length(x) == length(y)

Conclusions

In summary: if you use zip always ask yourself a question if you assume that the passed
iterators have the same length. If you do – make sure to check it in your code explicitly.