By: Logan Kilpatrick
Re-posted from: https://towardsdatascience.com/why-julia-2-0-isnt-coming-anytime-soon-and-why-that-is-a-good-thing-641ae3d2a177?source=rss-2c8aac9051d3------2
If you have been programming for a while, there’s a good chance you have lived through breaking changes in ecosystems which may have left…
Re-posted from: https://bkamins.github.io/julialang/2022/09/09/ssc2022.html
During SSC2022 Conference together with Rajith Vidanaarachchi and
Przemysław Szufel I am going to give
An introduction to agent-based simulations in the Julia language
tutorial.
During the tutorial tutorial agent-based models implemented in Julia will be
discussed. Today I thought to post an example of a simple agent-based model
(not covered during SSC2022) that shows one of the challenges of designing
agent-based models.
The post was written under Julia 1.8.0, DataFrames.jl 1.3.5, and Plots.jl 1.32.0.
I want to create a simulation of genetic drift in Julia.
The idea of the model is as follows. On a rectangular grid we initially randomly
put agents of two types (I will use true and false to denote them). In one
step we pick a random agent and change its type to the type of one of its
neighbors. We will want to check, using simulation, if in the long run all
agents end up having the same color and how long it takes, on the average,
to reach this state.
A crucial design decision when modeling this problem is how we define neighbors
of an agent. Some classical definitions are Moore or von Neumann
neighborhoods on a torus. However, in the experiment today I assume
that if I have an agent in location (x,y) it has four neighbors on diagonals:
(x+1, y+1), (x-1, y+1), (x+1, y-1), (x-1, y-1) (we assume the simulation
is run on a torus, so left and right edge, and top and bottom edges of the grid
are glued together). I want to see what consequence this definition of the
neighborhood has on the results of the simulation.
Let us start with an implementation on a grid of size 10 times 10 and plot the
dynamics of the system. Here is the code running 30,000 steps of the simulation
and want to see if we end up with all cells having the same type:
using Random
using Plots
Random.seed!(1234)
m = rand(Bool, 10, 10)
anim = @animate for tick in 1:30000
x, y = rand(axes(m, 1)), rand(axes(m, 2))
nx = mod1(x + rand((1, -1)), size(m, 1))
ny = mod1(y + rand((1, -1)), size(m, 2))
m[x, y] = m[nx, ny]
heatmap(m, legend=false, size=(400, 400), title=tick)
end every 100
gif(anim, "anim_fps15.gif", fps=15)
The code saves the result as an animated GIF file that you can see below:
Surprisingly, we reach a steady state that consists of 50% of true and 50%
of false cells. This teaches us that one should be careful when designing
ABM models, as small changes in the code can lead to significant changes it
the behavior of the model. The point is that our neighborhood consisting of
four diagonal neighbors defines two disjoint sets of cells on a square having even
size. In von Neumann neighborhood, which also consists of four negihbors
(but top-down and left-right) we would not have such a problem.
Let us investigate our flawed model for different sizes of the grid.
We want to check how often different steady states are reached for different
sizes of the grid and how long it takes to reach these steady states.
Before we can do it we first need to think of some rule that detects that
the simulation has reached a steady state and at the same time is cheap to check.
There are many ways to do it. Let me propose the following reasoning here. If we
have a square toroidal grid whose side is d we have 2d^2 pairs of nodes that
are neighbors. To see this note that we have d^2 agents and each node has 4
neighbors; in other words the grid is a 4-regular graph on d^2 nodes, so it
has 2d^2 edges.
Notice that if we in a continuous sequence visit every edge of this graph at
least once and we always see the same color on both ends of the edge we have
reached a steady state. Additionally, each edge of the graph is visited with
the same probability in each tick. Now assume that we run this process for t
ticks. The probability that we have never picked some concrete edge during this
time at is (1-1/(2d^2))^t. Since we have 2d^2 such edges, by union bound,
the probability that we have not picked at least one of the edges is less than
2d^2(1-1/(2d^2))^t. We want t to be large enough so that this probability is
small. Assume that we want this probability to be less than 1e-6. We thus
need to find a minimum t that satisfies 2d^2(1-1/(2d^2))^t<1e-6 as a
function of d. Using the approximation log(1+x) ≈ x for x close to 0 we
get that t > 2d^2*log(500_000d^2). For example for d=10 we can compute that
t must be at least 3546. In other words if for 3546 steps we do not make
any swaps of agent’s state we are very likely to have reached a steady state.
Using this derivation let us put our simulation in a function:
function runsim(d)
m = rand(Bool, d, d)
tick = 0
lastchange = 0
t = 2d^2*log(500_000d^2)
while lastchange < t
tick += 1
x, y = rand(axes(m, 1)), rand(axes(m, 2))
nx = mod1(x + rand((1, -1)), size(m, 1))
ny = mod1(y + rand((1, -1)), size(m, 2))
old = m[x, y]
new = m[nx, ny]
if old == new
lastchange += 1
else
lastchange = 0
end
m[x, y] = new
end
class=sum(m)
if iseven(d)
@assert allequal(m[i, j] for i in 1:d, j in 1:d if iseven(i+j))
@assert allequal(m[i, j] for i in 1:d, j in 1:d if isodd(i+j))
else
@assert allequal(m)
end
return (;d, tick, class)
end
Notice, that I have additionally added @assert checks that make sure
that indeed we have reached a steady state. However, we make this check only
once as it is expensive.
In the output the class value is number of true cells in the grid.
Let us now run this simulation for d ranging from 4 to 11 and 2048 times
for each size of the grid and keep the results in a data frame:
julia> using DataFrames
julia> using Statistics
julia> df = DataFrame([runsim(d) for d in 4:11 for _ in 1:2048])
16384×3 DataFrame
Row │ d tick class
│ Int64 Int64 Int64
───────┼─────────────────────
1 │ 4 595 8
2 │ 4 556 8
3 │ 4 570 8
⋮ │ ⋮ ⋮ ⋮
16382 │ 11 16112 0
16383 │ 11 17230 121
16384 │ 11 20914 0
16378 rows omitted
First, we want to analyze if indeed for even grid size we have three different
possible outcomes and for odd sized grid we have two possible steady states:
julia> combine(groupby(df, [:d, :class], sort=true), nrow, :tick .=> [mean, std])
20×5 DataFrame
Row │ d class nrow tick_mean tick_std
│ Int64 Int64 Int64 Float64 Float64
─────┼────────────────────────────────────────────
1 │ 4 0 504 604.919 68.1656
2 │ 4 8 1017 605.345 65.5917
3 │ 4 16 527 607.355 69.8784
4 │ 5 0 1015 1309.11 388.342
5 │ 5 25 1033 1341.45 415.868
6 │ 6 0 528 1853.26 425.983
7 │ 6 18 1038 1844.34 409.991
8 │ 6 36 482 1879.57 458.587
9 │ 7 0 1002 4077.81 2059.13
10 │ 7 49 1046 3901.26 1777.57
11 │ 8 0 521 4739.05 1562.39
12 │ 8 32 1027 4620.13 1481.39
13 │ 8 64 500 4717.03 1626.72
14 │ 9 0 1025 9825.73 5347.35
15 │ 9 81 1023 10069.5 5665.56
16 │ 10 0 491 10331.5 4139.39
17 │ 10 50 1075 10233.4 4008.2
18 │ 10 100 482 10469.8 4337.41
19 │ 11 0 1035 21127.0 12885.3
20 │ 11 121 1013 21374.3 13331.7
Indeed this is a case. Moreover we note that:
d has approximately the sameFinally let us calculate the average time to reach steady state by d and plot
it:
julia> res = combine(groupby(df, :d, sort=true), :tick => mean)
8×2 DataFrame
Row │ d tick_mean
│ Int64 Float64
─────┼──────────────────
1 │ 4 605.757
2 │ 5 1325.42
3 │ 6 1854.93
4 │ 7 3987.64
5 │ 8 4674.04
6 │ 9 9947.51
7 │ 10 10312.6
8 │ 11 21249.3
julia> scatter(res.d, res.tick_mean, legend=nothing, xlab="d", ylab="mean ticks")
The resulting plot is:
What is interesting is that there is much bigger increase of the time when we
move from even to odd d than when moving from odd to even d. What is the
reason for this? The reason is that for odd d the agents’ graph is connected.
There is a single component having d^2 elements. While for even d we have
two connected components, each having d^2/2 size. We have just learned that it
is much faster to independently reach a steady state for component of size
d^2/2 twice than to reach a steady state of a single component of size d^2.
I hope you enjoyed the post both from ABM and Julia language perspectives. I
tried to avoid going into too much detail, both from mathematical and
implementation side of the model. If you have any questions please contact me
on Julia Slack or Julia Discourse.
Re-posted from: https://bkamins.github.io/julialang/2022/09/02/speed.html
In my post from last week I tried to convince novice Julia users
to focus on writing simple code, and only when they learn more to start
writing fully generic code.
This week I continue my simple is best miniseries of posts. This time
I want to argue that you do not need to try writing code that
is maximally fast. Usually it is enough that your code is just fast
(fortunately in Julia, if you do not do some serious mistake your code
is not likely to be slow). Instead, I recommend that you try your code
to be simple and follow common patterns. In the long run it will be
much easier to read and update such code.
For the post I use DataFrames.jl as an example. It was tested under Julia 1.8.0,
DataFrames.jl 1.3.4, ShiftedArrays.jl 1.0.0, and BenchmarkTools.jl 1.3.1.
Recently one of the users asked the following question. Given a data frame
with column col add 5 columns to it with lags of column col ranging from
1 to 5.
Here is an example how you can do this task:
julia> using DataFrames
julia> using ShiftedArrays
julia> using BenchmarkTools
julia> df = DataFrame(a=1:10)
10×1 DataFrame
Row │ a
│ Int64
─────┼───────
1 │ 1
2 │ 2
3 │ 3
4 │ 4
5 │ 5
6 │ 6
7 │ 7
8 │ 8
9 │ 9
10 │ 10
julia> function add_lags(df)
out_df = copy(df)
for i in 1:5
out_df[:, "a$i"] = lag(df.a, i)
end
return out_df
end
add_lags (generic function with 1 method)
julia> add_lags(df)
10×6 DataFrame
Row │ a a1 a2 a3 a4 a5
│ Int64 Int64? Int64? Int64? Int64? Int64?
─────┼────────────────────────────────────────────────────
1 │ 1 missing missing missing missing missing
2 │ 2 1 missing missing missing missing
3 │ 3 2 1 missing missing missing
4 │ 4 3 2 1 missing missing
5 │ 5 4 3 2 1 missing
6 │ 6 5 4 3 2 1
7 │ 7 6 5 4 3 2
8 │ 8 7 6 5 4 3
9 │ 9 8 7 6 5 4
10 │ 10 9 8 7 6 5
julia> @btime add_lags($df);
2.900 μs (50 allocations: 3.26 KiB)
I have additionally added timing of this operation for a reference.
The add_lags function is meant to be an example of “fast code”
(it could still be made faster, or I could have chosen a more complex example
but I wanted to pick something that is easy to understand).
The functions combine, select[!], transform[!], and subset[!] support
the operation specification syntax
and are called the high-level API in DataFrames.jl.
A natural code for adding lags using the transform function is:
julia> transform(df, ["a" => (x -> lag(x, i)) => "a$i" for i in 1:5])
10×6 DataFrame
Row │ a a1 a2 a3 a4 a5
│ Int64 Int64? Int64? Int64? Int64? Int64?
─────┼────────────────────────────────────────────────────
1 │ 1 missing missing missing missing missing
2 │ 2 1 missing missing missing missing
3 │ 3 2 1 missing missing missing
4 │ 4 3 2 1 missing missing
5 │ 5 4 3 2 1 missing
6 │ 6 5 4 3 2 1
7 │ 7 6 5 4 3 2
8 │ 8 7 6 5 4 3
9 │ 9 8 7 6 5 4
10 │ 10 9 8 7 6 5
The problem with this solution that users often raise is that it is slower than
the previous one:
julia> @btime transform($df, ["a" => (x -> lag(x, i)) => "a$i" for i in 1:5]);
61.600 μs (523 allocations: 27.46 KiB)
There are two questions that one naturally asks. The first is what is the reason
of this timing difference and the second is do we care.
So first let me explain the reason for the slowdown. The transform function is
flexible and allows for many different operations. Therefore it does,
apart from doing core transformations, a lot of extra pre and post processing
that is needed to provide this flexibility.
Now let us try to answer the question if we care. First check a bigger data frame:
julia> df2 = DataFrame(a=1:10^8);
julia> @btime add_lags($df2);
1.354 s (62 allocations: 4.94 GiB)
julia> @btime transform($df2, ["a" => (x -> lag(x, i)) => "a$i" for i in 1:5]);
1.318 s (539 allocations: 4.94 GiB)
The timings are roughly the same.
As you can see – the pre and post processing time in transform does not grow
with size of data. Therefore we get the following conclusions:
transform might end up being a significantYou might think that the code:
["a" => (x -> lag(x, i)) => "a$i" for i in 1:5]
is not so easy to read. This might be true for a newcomer, as this is something
new you need to learn. However, my experience is that after some practice with
DataFrames.jl operation specification language it becomes quite easy to write
and read. You just need to embrace the pattern:
[input column] => [transformation function] => [output column name]
and the rest is a natural consequence (one of the benefits of this pattern
is that it is clear visually as => nicely separates its parts).
The benefit of the high-level API is that most likely the lagging operation we
have considered needs to be done in groups. That is, you have a grouping
variable, call it id, and you want to perform lagging per-id value
(e.g. you have different stocks and their quotations for different days
that you want to lag).
With transform this is a breeze. You just need to add groupby to your code:
julia> df3 = DataFrame(id=repeat(1:2, inner=8), a=1:16)
16×2 DataFrame
Row │ id a
│ Int64 Int64
─────┼──────────────
1 │ 1 1
2 │ 1 2
3 │ 1 3
4 │ 1 4
5 │ 1 5
6 │ 1 6
7 │ 1 7
8 │ 1 8
9 │ 2 9
10 │ 2 10
11 │ 2 11
12 │ 2 12
13 │ 2 13
14 │ 2 14
15 │ 2 15
16 │ 2 16
julia> transform(groupby(df3, :id),
["a" => (x -> lag(x, i)) => "a$i" for i in 1:5])
16×7 DataFrame
Row │ id a a1 a2 a3 a4 a5
│ Int64 Int64 Int64? Int64? Int64? Int64? Int64?
─────┼───────────────────────────────────────────────────────────
1 │ 1 1 missing missing missing missing missing
2 │ 1 2 1 missing missing missing missing
3 │ 1 3 2 1 missing missing missing
4 │ 1 4 3 2 1 missing missing
5 │ 1 5 4 3 2 1 missing
6 │ 1 6 5 4 3 2 1
7 │ 1 7 6 5 4 3 2
8 │ 1 8 7 6 5 4 3
9 │ 2 9 missing missing missing missing missing
10 │ 2 10 9 missing missing missing missing
11 │ 2 11 10 9 missing missing missing
12 │ 2 12 11 10 9 missing missing
13 │ 2 13 12 11 10 9 missing
14 │ 2 14 13 12 11 10 9
15 │ 2 15 14 13 12 11 10
16 │ 2 16 15 14 13 12 11
The add_lags function does not obviously allow you to do what you want.
After some thinking (and if you know DataFrames.jl well enough) you come to the
conclusion that the following will work:
julia> transform(groupby(df3, :id), add_lags)
16×7 DataFrame
Row │ id a a1 a2 a3 a4 a5
│ Int64 Int64 Int64? Int64? Int64? Int64? Int64?
─────┼───────────────────────────────────────────────────────────
1 │ 1 1 missing missing missing missing missing
2 │ 1 2 1 missing missing missing missing
3 │ 1 3 2 1 missing missing missing
4 │ 1 4 3 2 1 missing missing
5 │ 1 5 4 3 2 1 missing
6 │ 1 6 5 4 3 2 1
7 │ 1 7 6 5 4 3 2
8 │ 1 8 7 6 5 4 3
9 │ 2 9 missing missing missing missing missing
10 │ 2 10 9 missing missing missing missing
11 │ 2 11 10 9 missing missing missing
12 │ 2 12 11 10 9 missing missing
13 │ 2 13 12 11 10 9 missing
14 │ 2 14 13 12 11 10 9
15 │ 2 15 14 13 12 11 10
16 │ 2 16 15 14 13 12 11
so things are still easy to code (but need transform).
Let us compare the performance of both operations on a larger data frame:
julia> df4 = DataFrame(id=repeat(1:10, inner=10^7), a=1:10^8);
julia> @btime transform(groupby($df4, :id),
["a" => (x -> lag(x, i)) => "a$i" for i in 1:5]);
10.124 s (1731 allocations: 19.35 GiB)
julia> @btime transform(groupby($df4, :id), add_lags);
11.747 s (1552 allocations: 24.59 GiB)
We see that using add_lags in this case is slightly slower. Of course we could
write a custom add_lags function that would work by groups, even without using
groupby but this would be not so easy (I recommend you try it if you are not
convinced).
Here we see the benefit of general design of transform: it handles grouped
data in the same way as it handles a data frame.
Let me summarize my thinking about code performance in Julia in general,
and in DataFrames.jl in particular:
transform in our case) do will be negligible anyway (this