# Simulating billiards in Julia

<p class="pm-node nj-authors">George Datseris, a[link]</p>

# Author

[George Datseris](https://github.com/Datseris), contributor of [JuliaDynamics](https://github.com/JuliaDynamics) and [JuliaMusic](https://github.com/JuliaDynamics)

# This tutorial

In this tutorial we will simulate a dynamical billiard. This is a very simple system where a point particle is propagating inside a domain with constant speed. When encountering a boundary, the particle undergoes specular reflection.

# Features of Julia highlighted

* Intuitive coding of the problem
* High-end Performance without doing anything
* Simple metaprogramming for performance gain
* Extendability with least possible effort
* Extandability applies to *all* aspects of the system (Multiple Dispatch)

# The system and algorithm

* Particles that propagate inside the billiard
* Obstacles that compose the billiard
* Processes that perform the propagation

# The Algorithm

* Calculate the collisions of `p` with all obstacles in the billiard
* Find the collision that happens first in time, and the obstacle it corresponds to
* Propagate the particle to this collision point
* Do specular reflection at the boundary of the obstacle to be collided with
* Rinse and repeat!

# The basic structs we will need

All vectors we'll use are two dimensional. We can therefore take advantage of `StaticArrays`. For a shortcut, I'll define `SV` to be a 2D vector of Floats (all the vectors we will use will be of this type)

```julia id=8e2af7db-bc64-40d0-a99d-4619558f7c51
using StaticArrays
const SV = SVector{2,Float64}
```

We now define the particle `struct`. For extendability, it is best to define it based on an abstract type.

```julia id=923873c8-b3da-495b-a252-8326d6e34bd1
abstract type AbstractParticle end

mutable struct Particle <: AbstractParticle
    pos::SV
    vel::SV
end
Particle(x0, y0, φ0) = Particle(SV(x0, y0), SV(cos(φ0), sin(φ0)))
```

Then, we define the obstacle `struct`s. Again, to allow extendability, abstraction and generalization, it is best to define all obstacles as a subtype of an abstract type.

```julia id=769c00aa-ab1b-4129-baf4-6010e5fc1a84
abstract type Obstacle end

struct Wall <: Obstacle
    sp::SV
    ep::SV
    normal::SV
end

struct Disk <: Obstacle
    c::SV
    r::Float64
end
```

A "billiard" in the following will simply be a `Tuple` of `Obstacle` subtypes. Specifically, the type of the billiard is `NTuple{N, Obstacle} where {N}`, so it is convenient to define:

```julia id=66e1e46d-dab7-412c-a250-c62de6898505
const Billiard = NTuple{N, Obstacle} where N
```

# Collisions

As said in the introduction, the most important part is a function that finds collisions between particles and obstacles. Here is how this looks for collision of a standard particle with a wall:

```julia id=64ef2a21-ae5f-415f-b193-8d3d2f26070f
using LinearAlgebra: dot, normalize

"""
    collision(p::AbstractParticle, o::Obstacle) → t, cp
Find the collision (if any) between given particle and obstacle.
Return the time until collision and the estimated collision point `cp`.
"""
@inline function collision(p::Particle, w::Wall)
    n = normalvec(w, p.pos)
    denom = dot(p.vel, n)
    if denom ≥ 0.0
        return nocollision()
    else
        t = dot(w.sp - p.pos, n)/denom
        return t, p.pos + t * p.vel
    end
end

normalvec(w::Wall, pos) = w.normal
```

and here is the same function but for collisions with disks instead:

```julia id=44dc5b78-8e8a-4db4-bd83-735abe8c699e
@inline function collision(p::Particle, d::Disk)
    dotp = dot(p.vel, normalvec(d, p.pos))
    dotp ≥ 0.0 && return nocollision()

    dc = p.pos - d.c
    B = dot(p.vel, dc)           #pointing towards circle center: B < 0
    C = dot(dc, dc) - d.r*d.r    #being outside of circle: C > 0
    Δ = B*B - C

    Δ ≤ 0.0 && return nocollision()
    sqrtD = sqrt(Δ)
    # Closest point:
    t = -B - sqrtD
    return t, p.pos + t * p.vel
end

normalvec(d::Disk, pos) = normalize(pos - d.c)
```

You can see that there are cases where collisions are not possible *or* they happen backwards in time. By convention, this is the value we return then:

```julia id=8a50c18c-1ae3-48c2-bd44-d212807826c0
@inline nocollision() = (Inf, SV(0.0, 0.0))
```

`next_collision` is a useful function that finds the "true" next collision. It simply loops over the obstacles in a billiard. It simply checks which obstacle has the least collision time:

```julia id=52fabede-b7f7-47f3-8582-1fa621568f6f
function next_collision(p::AbstractParticle, bd)
    j, ct, cp = 0, Inf, SV(0.0, 0.0)
    for i in eachindex(bd)
        t, c = collision(p, bd[i])
        if t < ct
            j = i
            ct = t
            cp = c
        end
    end
    return j, ct, cp
end
```

# Evolving a particle in a billiard

We need a simple function to propagate a particle to the found collision point. We will also a give the amount of time the propagation "should" take, for extendability. For standard particles, where the velocity vector is constant while travelling, this does not matter.

```julia id=6f1e702d-f90b-4307-a1de-2de8ac296e90
propagate!(p::Particle, pos, t) = (p.pos = pos)
```

This is an *in-place* function (notice `!` at the end)

```julia id=9e16989a-cc81-4bf7-884c-7be652ba07c4
function resolvecollision!(p::AbstractParticle, o::Obstacle)
    n = normalvec(o, p.pos)
    p.vel = p.vel - 2*dot(n, p.vel)*n
end
```

We are now ready to wrap things up. Let's define a function that takes a particle and evolves in a billiard (tuple of obstacles) and returns the timeseries of the positions of the particle.

```julia id=bf4c6776-58c1-485e-ac36-e52d4c35d037
"""
    bounce!(p, bd)
Evolve the particle for one collision (in-place).
"""
@inline function bounce!(p::AbstractParticle, bd)
    i::Int, tmin::Float64, cp::SV = next_collision(p, bd)
    if tmin != Inf
        propagate!(p, cp, tmin)
        resolvecollision!(p, bd[i])
    end
    return i, tmin, p.pos, p.vel
end
```

Then, we can use it inside a bigger function that calls `bounce!` until a specified amount of time:

```julia id=8fcfde6b-f9d7-481c-9b83-c6d9748c0e33
"""
    timeseries!(p::AbstractParticle, bd, n) -> xt, yt, t
Evolve the particle in the billiard `bd` for `n` collisions
and return the position timeseries `xt, yt` along with time vector `t`.
"""
function timeseries!(p::AbstractParticle, bd, n::Int)

    t = [0.0]; xt = [p.pos[1]]; yt = [p.pos[2]]; c = 0

    while c < n
        prevpos = p.pos; prevvel = p.vel
        i, ct = bounce!(p, bd)
        xs, ys = extrapolate(p, prevpos, prevvel, ct)
        push!(t, ct)
        append!(xt, xs)
        append!(yt, ys)
        c += 1
    end

    return xt, yt, t
end
```

`extrapolate` simply creates the position timeseries in between two collisions. For a standard particle there is no "extrapolation" needed, one just uses the final position:

```julia id=20ee2a53-e9b2-4718-ace0-8ac6ccb71c29
extrapolate(p::Particle, prevpos, prevvel, ct) = p.pos
```

Why does this `extrapolate` function exist? See below when we extend our code for magnetic particles!

# Running the code

let's put this to the test now! We'll create the famous Sinai billiard

```julia id=945798fe-6038-4c34-8178-7be7d4048f9a
x, y, r = 1.0, 1.0, 0.3
sp = [0.0,y]; ep = [0.0, 0.0]; n = [x,0.0]
leftw = Wall(sp, ep, n)
sp = [x,0.0]; ep = [x, y]; n = [-x,0.0]
rightw = Wall(sp, ep, n)
sp = [x,y]; ep = [0.0, y]; n = [0.0,-y]
topw = Wall(sp, ep, n)
sp = [0.0,0.0]; ep = [x, 0.0]; n = [0.0,y]
botw = Wall(sp, ep, n)
disk = Disk([x/2, y/2], r)

bd = (botw, rightw, topw, leftw, disk)
bd isa Billiard
```

and also initialize a particle

```julia id=34b9bf2f-e00d-4b8a-bd48-222a5e833e36
p = Particle(0.1, 0.1, 2π*rand())
```

and evolve it

```julia id=4d36cc2e-9dd1-483e-9ee8-75d0905d4a2a
xt, yt, t = timeseries!(p, bd, 10)
```

# Plotting

Let's define some simple methods for plotting and plot the result!

```julia id=25d863da-8535-4a19-8f6e-31ca470e5eee
using PyPlot
import PyPlot: plot

const EDGECOLOR = (0,0.6,0)
function plot(d::Disk)
    facecolor = (EDGECOLOR..., 0.5)
    circle1 = PyPlot.plt[:Circle](d.c, d.r;
        edgecolor = EDGECOLOR, facecolor = facecolor, lw = 2.0)
    PyPlot.gca()[:add_artist](circle1)
end
function plot(w::Wall)
    PyPlot.plot([w.sp[1],w.ep[1]],[w.sp[2],w.ep[2]];
        color=EDGECOLOR, lw = 2.0)
end
function plot(bd::Billiard)
    for o ∈ bd; plot(o); end
    gca()[:set_aspect]("equal")
end

figure(); plot(bd)
gcf()
```

Awesome! Now let's see it with the orbit as well

```julia id=eb21af35-6b91-4601-82b2-ac95369a1d36
figure(); plot(bd)
xt, yt, t = timeseries!(p, bd, 10)
plot(xt, yt)
gcf()
```

Plot 10 orbits!

```julia id=e9219790-56fd-49ae-8a30-8a9e998f69e4
figure(); plot(bd)
p = Particle(0.1, 0.5, 2π*rand())
for j in 1:10
    xt, yt = timeseries!(p, bd, 20)
    plot(xt, yt, alpha = 0.5)
end
gcf()
```

# Showcase 1: Performance & Metaprogramming

```julia id=eb750630-d440-4ae8-8492-0ed9930e914f
using BenchmarkTools

p = Particle(0.1, 0.1, 2π*rand())
@btime bounce!($p, $bd)
```

It is already very fast to propagate a particle for one collision, however there are some allocations (even though the function is in place).

```julia id=82707cc9-243d-4a69-abd2-b49cc8f4fc19
@generated function next_collision(p::AbstractParticle, bd::Billiard)
    L = length(bd.types) # notice that bd stands for the TYPE of bd here!
    out = :(ind = 0; tmin = Inf; cp = SV(0.0, 0.0))
    for j=1:L
        push!(out.args, quote
                            let x = bd[$j]
                                tcol, pcol = collision(p, x)
                                # Set minimum time:
                                if tcol < tmin
                                  tmin = tcol
                                  ind = $j
                                  cp = pcol
                                end
                            end
                        end
                        )
    end
    push!(out.args, :(return ind, tmin, cp))
    return out
end


@btime bounce!($p, $bd)
```

This number is **insane**!!! Notice that this code is billiard agnostic! You could pass any tuple of obstacles and it would still be as performant!!! The time of `bounce!` scales linearly with the number of obstacles in the billiard.

# Showcase 2: Extendability

Let's say we want to add one more obstacle to this "billiard package" we are making. Do you we have to re-write *everything* for it? Nope! In the end we only need to extend two methods ! **Only two**!

```julia id=fc4427ba-0ce5-4940-93a2-0b5796532d09
struct Ellipse <: Obstacle
    c::SV
    a::Float64
    b::Float64
end
```

The methods we need to extend are only these:

```julia id=ae70352a-6d75-49f7-a8cc-793b9f43008c
normalvec
collision
```

Yes!!! Only two! So let's get to it! `normalvec` is pretty easy:

```julia id=9d1b09a9-b0f3-4ab7-a183-979565ba7871
function normalvec(e::Ellipse, pos)
    x₀, y₀ = pos
    h, k = e.c
    return normalize(SV((x₀-h)/(e.a*e.a), (y₀-k)/(e.b*e.b)))
end

using LinearAlgebra: norm
function collision(p::Particle, e::Ellipse)
    dotp = dot(p.vel, normalvec(e, p.pos))
    dotp ≥ 0.0 && return nocollision()

    a = e.a; b = e.b
    pc = p.pos - e.c
    μ = p.vel[2]/p.vel[1]
    ψ = pc[2] - μ*pc[1]

    denomin = a*a*μ*μ + b*b
    Δ² = denomin - ψ*ψ
    Δ² ≤ 0 && return nocollision()
    Δ = sqrt(Δ²); f1 = -a*a*μ*ψ; f2 = b*b*ψ # just factors
    I1 = SV(f1 + a*b*Δ, f2 + a*b*μ*Δ)/denomin
    I2 = SV(f1 - a*b*Δ, f2 - a*b*μ*Δ)/denomin

    d1 = norm(pc - I1); d2 = norm(pc - I2)
    return d1 < d2 ? (d1, I1 + e.c) : (d2, I2 + e.c)
end
```

Alright so now let's create a billiard with both an ellipse and a disk, for the fun of it

```julia id=5c42e374-2ff1-4603-804c-03b349df69de
el = Ellipse([0.4, 0.2 ], 0.3, 0.1)
di = Disk([0.6, 0.7], 0.25)

bd2 = Billiard((bd[1:4]..., el, di))
```

and plot it

```julia id=4cfeb44a-3686-4268-9464-273ada10f857
function plot(e::Ellipse)
    facecolor = (EDGECOLOR..., 0.5)
    ellipse = PyPlot.matplotlib[:patches][:Ellipse](e.c, 2e.a, 2e.b;
        edgecolor = EDGECOLOR, facecolor = facecolor, lw = 2.0)
    PyPlot.gca()[:add_artist](ellipse)
end

figure(); plot(bd2)
gcf()
```

We are now ready to evolve a particle in this brand new billiard:

```julia id=c0ef7b56-0db3-42c7-be67-66873ccc38f5
p = Particle(0.1, 0.1, 2π*rand())
xt, yt, t = timeseries!(p, bd2, 20)
figure(); plot(bd2)
plot(xt, yt)
gcf()
```

plot a bunch more!

```julia id=70a1772c-ba02-465e-b1a2-e3dd90869363
figure(); plot(bd2)
for j in 1:10
    p = Particle(0.1, 0.1, 2π*rand())
    xt, yt = timeseries!(p, bd2, 20)
    plot(xt, yt, alpha = 0.5)
end
gcf()
```

# Showcase 3: Extendability, again.

Alright, so it turned out to be almost trivial to add an extra obstacle to our code. But what about an extra particle?

```julia id=e9d06f18-94ef-4cb0-aa39-384afed4cde1
collision # for each obstacle we want to support
propagate!
extrapolate
```

and yeap, that's it. It may be hard to believe that it only takes so little, but it's true!!!

## The type

```julia id=0bc7d0ce-affd-4e23-92b7-2378f67a427c
mutable struct MagneticParticle <: AbstractParticle
    pos::SV
    vel::SV
    ω::Float64
end
MagneticParticle(x0, y0, φ0, ω) = MagneticParticle(SV(x0, y0), SV(cos(φ0), sin(φ0)), ω)
```

This particle moves in circles with angular velocity `ω`.

## Extending collision

To extend `collision`, we simply have to find intersections of circle-line and circle-circle, for collisions with `Wall` and `Disk`. I won't go into details of how to do this, and instead I'll copy-paste functions from `DynamicalBilliards`. The versions in `DynamicalBilliards` also have a lot of comments that explain what is going on.

```julia id=f9d4519f-e442-4d01-86f5-b85af86ab24d
function collision(p::MagneticParticle, w::Wall)
    ω = p.ω
    pc, pr = cyclotron(p)
    P0 = p.pos
    P2P1 = w.ep - w.sp
    P1P3 = w.sp - pc

    a = dot(P2P1, P2P1)
    b = 2*dot(P2P1, P1P3)
    c = dot(P1P3, P1P3) - pr*pr
    Δ = b^2 -4*a*c

    Δ ≤ 0.0 && return nocollision()
    u1 = (-b - sqrt(Δ))/2a
    u2 = (-b + sqrt(Δ))/2a
    cond1 = 0.0 ≤ u1 ≤ 1.0
    cond2 = 0.0 ≤ u2 ≤ 1.0
    θ, I = nocollision()
    if cond1 || cond2
        dw = w.ep - w.sp
        for (u, cond) in ((u1, cond1), (u2, cond2))
            Y =  w.sp + u*dw
            if cond
                φ = realangle(p, w, Y)
                φ < θ && (θ = φ; I = Y)
            end
        end
    end
    return θ*pr, I
end
```

and here is the collision with a disk:

```julia id=b7dd83f0-83e6-431e-9949-73ccb692be6c
function collision(p::MagneticParticle, o::Disk)
    ω = p.ω
    pc, rc = cyclotron(p)
    p1 = o.c
    r1 = o.r
    d = norm(p1-pc)
    if (d >= rc + r1) || (d <= abs(rc-r1))
        return nocollision()
    end

    a = (rc^2 - r1^2 + d^2)/2d
    h = sqrt(rc^2 - a^2)

    I1 = SV(
        pc[1] + a*(p1[1] - pc[1])/d + h*(p1[2] - pc[2])/d,
        pc[2] + a*(p1[2] - pc[2])/d - h*(p1[1] - pc[1])/d
    )
    I2 = SV(
        pc[1] + a*(p1[1] - pc[1])/d - h*(p1[2] - pc[2])/d,
        pc[2] + a*(p1[2] - pc[2])/d + h*(p1[1] - pc[1])/d
    )

    θ1 = realangle(p, o, I1)
    θ2 = realangle(p, o, I2)

    return θ1 < θ2 ? (θ1*rc, I1) : (θ2*rc, I2)
end
```

The functions `cyclotron` and `realangle` are helper functions. The first one finds the center and radius of the cyclotron traced by the particle.

```julia id=dfdff7c7-9611-4ae8-b9fd-f3f1828a1de4
cyclotron(p) = (p.pos - (1/p.ω)*SV(p.vel[2], -p.vel[1]), abs(1/p.ω))
```

`realangle` has a simple purpose: the intersections of a circle with any obstacle are always 2. But which one happens first, from a temporal perspective? `realangle` gives the correct angle until the collision point, in forward time.

```julia id=a79ae0f2-ed5d-4980-b587-b4572dbc0c7a
function realangle(p::MagneticParticle, o::Obstacle, i)

    pc, pr = cyclotron(p); ω = p.ω
    P0 = p.pos
    PC = pc - P0
    d2 = dot(i-P0,i-P0)

    if d2 ≤ 1e-8
        dotp = dot(p.vel, normalvec(o,  p.pos))
        dotp ≥ 0 && return Inf
    end

    d2r = (d2/(2pr^2))
    d2r > 2 && (d2r = 2.0)
    θprime = acos(1.0 - d2r)

    PI = i - P0
    side = (PI[1]*PC[2] - PI[2]*PC[1])*ω

    side < 0 && (θprime = 2π-θprime)
    return θprime
end
```

The complexity of the functions `collision` and `realangle` exists solely due to the geometry of intersections between circles. What we want to point out is how **few** methods we have to extend. How easy is defining these new methods is not relevant, blame math and physics for that! So don't be taken aback because these functions are "long"!

## Propagation & extrapolation

`propagate!` for a `MagneticParticle` must evolve it in an arc of a circle, so as you can see we have to change the velocity vector!

```julia id=87d1a870-d414-4f41-9bf6-661f6106d1b9
function propagate!(p::MagneticParticle, pos, t)
    φ0 = atan(p.vel[2], p.vel[1])
    p.pos = pos
    p.vel = SV(cossin(p.ω*t + φ0))
    return
end

cossin(x) = ( (y, z) = sincos(x); (z, y) )
```

`extrapolate` should simply create the arc that connects the previous point with the current one

```julia id=35eef31e-b37e-44ef-8fdb-0d8f5bfbc671
function extrapolate(p::MagneticParticle, prevpos, prevvel, t)
    φ0 = atan(prevvel[2], prevvel[1])
    s0, c0 = sincos(φ0)
    x0 = prevpos[1]; y0 = prevpos[2]
    xt = [x0]; yt = [y0]; ω = p.ω
    tvec = 0.0:0.01:t
    for td in tvec
        s, c = sincos(p.ω*td + φ0)
        push!(xt, s/ω + x0 - s0/ω)
        push!(yt, -c/ω + y0 + c0/ω) #vx0 is cos(φ0)
    end
    return xt, yt
end
```

## Evolve the magnetic particle

```julia id=6bc1f289-1517-403b-9e1c-7ad725f9976f
p = MagneticParticle(0.1, 0.1, 2π*rand(), 3.0)
xt, yt, t = timeseries!(p, bd, 20)
figure(); plot(bd)
plot(xt, yt)
gcf()
```

plot a bunch of these!

```julia id=d4925d1a-dac9-4228-b469-7711c3b7d18e
figure(); plot(bd)
for j in 1:4
    p = MagneticParticle(0.1, 0.1, 2π*rand(), 2.0)
    xt, yt = timeseries!(p, bd, 20)
    plot(xt, yt, alpha = 0.5)
end
gcf()
```


<details id="com.nextjournal.article">
<summary>This notebook was exported from <a href="https://nextjournal.com/a/KQqarQPvsAfwYFPq1AadZ?change-id=CLPh3GwFARfcNn35Pf64H3">https://nextjournal.com/a/KQqarQPvsAfwYFPq1AadZ?change-id=CLPh3GwFARfcNn35Pf64H3</a></summary>

```edn nextjournal-metadata
{:article
 {:settings
  {:numbered? true,
   :sidebar? true,
   :authors? false,
   :centered? false,
   :subtitle? false},
  :nodes
  {"0bc7d0ce-affd-4e23-92b7-2378f67a427c"
   {:compute-ref #uuid "4f58c240-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 365,
    :id "0bc7d0ce-affd-4e23-92b7-2378f67a427c",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "20ee2a53-e9b2-4718-ace0-8ac6ccb71c29"
   {:compute-ref #uuid "3632ec50-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 272,
    :id "20ee2a53-e9b2-4718-ace0-8ac6ccb71c29",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "25d863da-8535-4a19-8f6e-31ca470e5eee"
   {:compute-ref #uuid "38a0b800-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 20660,
    :id "25d863da-8535-4a19-8f6e-31ca470e5eee",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "34b9bf2f-e00d-4b8a-bd48-222a5e833e36"
   {:compute-ref #uuid "36d2d5d0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 1059,
    :id "34b9bf2f-e00d-4b8a-bd48-222a5e833e36",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "35eef31e-b37e-44ef-8fdb-0d8f5bfbc671"
   {:compute-ref #uuid "508645c0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 198,
    :id "35eef31e-b37e-44ef-8fdb-0d8f5bfbc671",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "44dc5b78-8e8a-4db4-bd83-735abe8c699e"
   {:compute-ref #uuid "355ab4c0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 49,
    :id "44dc5b78-8e8a-4db4-bd83-735abe8c699e",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "4cfeb44a-3686-4268-9464-273ada10f857"
   {:compute-ref #uuid "4dc85e40-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 501,
    :id "4cfeb44a-3686-4268-9464-273ada10f857",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "4d36cc2e-9dd1-483e-9ee8-75d0905d4a2a"
   {:compute-ref #uuid "37749410-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 1966,
    :id "4d36cc2e-9dd1-483e-9ee8-75d0905d4a2a",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "52fabede-b7f7-47f3-8582-1fa621568f6f"
   {:compute-ref #uuid "357ede90-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 174,
    :id "52fabede-b7f7-47f3-8582-1fa621568f6f",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "5c42e374-2ff1-4603-804c-03b349df69de"
   {:compute-ref #uuid "4d5d6770-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 700,
    :id "5c42e374-2ff1-4603-804c-03b349df69de",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "64ef2a21-ae5f-415f-b193-8d3d2f26070f"
   {:compute-ref #uuid "338aafb0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 3040,
    :id "64ef2a21-ae5f-415f-b193-8d3d2f26070f",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "66e1e46d-dab7-412c-a250-c62de6898505"
   {:compute-ref #uuid "335017b0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 383,
    :id "66e1e46d-dab7-412c-a250-c62de6898505",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "6bc1f289-1517-403b-9e1c-7ad725f9976f"
   {:compute-ref #uuid "50a4a330-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 1171,
    :id "6bc1f289-1517-403b-9e1c-7ad725f9976f",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "6f1e702d-f90b-4307-a1de-2de8ac296e90"
   {:compute-ref #uuid "35999280-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 174,
    :id "6f1e702d-f90b-4307-a1de-2de8ac296e90",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "70a1772c-ba02-465e-b1a2-e3dd90869363"
   {:compute-ref #uuid "4ebc3420-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 944,
    :id "70a1772c-ba02-465e-b1a2-e3dd90869363",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "769c00aa-ab1b-4129-baf4-6010e5fc1a84"
   {:compute-ref #uuid "334e9110-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 9,
    :id "769c00aa-ab1b-4129-baf4-6010e5fc1a84",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "82707cc9-243d-4a69-abd2-b49cc8f4fc19"
   {:compute-ref #uuid "4a1911e0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 4536,
    :id "82707cc9-243d-4a69-abd2-b49cc8f4fc19",
    :kind "code",
    :output-log-lines {:stdout 1},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "87d1a870-d414-4f41-9bf6-661f6106d1b9"
   {:compute-ref #uuid "5049b1f0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 396,
    :id "87d1a870-d414-4f41-9bf6-661f6106d1b9",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "8a50c18c-1ae3-48c2-bd44-d212807826c0"
   {:compute-ref #uuid "356255e0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 186,
    :id "8a50c18c-1ae3-48c2-bd44-d212807826c0",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "8e2af7db-bc64-40d0-a99d-4619558f7c51"
   {:compute-ref #uuid "329ddeb0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 1107,
    :id "8e2af7db-bc64-40d0-a99d-4619558f7c51",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "8fcfde6b-f9d7-481c-9b83-c6d9748c0e33"
   {:compute-ref #uuid "362a87e0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 54,
    :id "8fcfde6b-f9d7-481c-9b83-c6d9748c0e33",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "923873c8-b3da-495b-a252-8326d6e34bd1"
   {:compute-ref #uuid "3346eff0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 49,
    :id "923873c8-b3da-495b-a252-8326d6e34bd1",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "945798fe-6038-4c34-8178-7be7d4048f9a"
   {:compute-ref #uuid "365c9460-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 774,
    :id "945798fe-6038-4c34-8178-7be7d4048f9a",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "9d1b09a9-b0f3-4ab7-a183-979565ba7871"
   {:compute-ref #uuid "4d3ac440-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 226,
    :id "9d1b09a9-b0f3-4ab7-a183-979565ba7871",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "9e16989a-cc81-4bf7-884c-7be652ba07c4"
   {:compute-ref #uuid "35b44670-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 274,
    :id "9e16989a-cc81-4bf7-884c-7be652ba07c4",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "a79ae0f2-ed5d-4980-b587-b4572dbc0c7a"
   {:compute-ref #uuid "500cd000-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 398,
    :id "a79ae0f2-ed5d-4980-b587-b4572dbc0c7a",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"
   {:environment
    [:environment
     {:article/nextjournal.id
      #uuid "02921106-cea7-469d-8e6d-fb4f112c7695",
      :change/nextjournal.id
      #uuid "5bc6060f-40b9-4f0e-ba3d-4e717ad8d6be",
      :node/id "2b507ea8-8459-4f12-bd3a-b084a13f0ea4"}],
    :id "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58",
    :kind "runtime",
    :language "julia",
    :type :nextjournal},
   "ae70352a-6d75-49f7-a8cc-793b9f43008c"
   {:compute-ref #uuid "4ccebc00-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 707,
    :id "ae70352a-6d75-49f7-a8cc-793b9f43008c",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "b7dd83f0-83e6-431e-9949-73ccb692be6c"
   {:compute-ref #uuid "4fae5c50-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 213,
    :id "b7dd83f0-83e6-431e-9949-73ccb692be6c",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "bf4c6776-58c1-485e-ac36-e52d4c35d037"
   {:compute-ref #uuid "35de3ca0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 499,
    :id "bf4c6776-58c1-485e-ac36-e52d4c35d037",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "c0ef7b56-0db3-42c7-be67-66873ccc38f5"
   {:compute-ref #uuid "4e14f7a0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 1095,
    :id "c0ef7b56-0db3-42c7-be67-66873ccc38f5",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "d4925d1a-dac9-4228-b469-7711c3b7d18e"
   {:compute-ref #uuid "c55c7080-d226-11e8-b7e3-d5a7297efcb7",
    :exec-duration 819,
    :id "d4925d1a-dac9-4228-b469-7711c3b7d18e",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "dfdff7c7-9611-4ae8-b9fd-f3f1828a1de4"
   {:compute-ref #uuid "4fcf03b0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 404,
    :id "dfdff7c7-9611-4ae8-b9fd-f3f1828a1de4",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "e9219790-56fd-49ae-8a30-8a9e998f69e4"
   {:compute-ref #uuid "457db7d0-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 1077,
    :id "e9219790-56fd-49ae-8a30-8a9e998f69e4",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "e9d06f18-94ef-4cb0-aa39-384afed4cde1"
   {:compute-ref #uuid "4f4c6630-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 81,
    :id "e9d06f18-94ef-4cb0-aa39-384afed4cde1",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "eb21af35-6b91-4601-82b2-ac95369a1d36"
   {:compute-ref #uuid "44f15650-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 919,
    :id "eb21af35-6b91-4601-82b2-ac95369a1d36",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "eb750630-d440-4ae8-8492-0ed9930e914f"
   {:compute-ref #uuid "46225c40-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 6649,
    :id "eb750630-d440-4ae8-8492-0ed9930e914f",
    :kind "code",
    :output-log-lines {:stdout 1},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "f9d4519f-e442-4d01-86f5-b85af86ab24d"
   {:compute-ref #uuid "4f90c230-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 193,
    :id "f9d4519f-e442-4d01-86f5-b85af86ab24d",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]},
   "fc4427ba-0ce5-4940-93a2-0b5796532d09"
   {:compute-ref #uuid "4ccd5c70-d225-11e8-b7e3-d5a7297efcb7",
    :exec-duration 9,
    :id "fc4427ba-0ce5-4940-93a2-0b5796532d09",
    :kind "code",
    :output-log-lines {},
    :runtime [:runtime "ac33edd0-9cb1-4d3d-a179-38d6e9b5eb58"]}},
  :nextjournal/id #uuid "0292140d-a5f1-496c-b3f3-3772026a6ad4",
  :article/change
  {:nextjournal/id #uuid "5bc9a3ea-ece9-496f-8a7c-4b35df6e92f6"}}}

```
</details>