3. Animations and Interaction with Traveling Salesperson Problem

In an attempt to try using a number of genetic algorithms packages I’ve been seeing, and also get more comfortable with WGLMakie and Bonito, I’ve decided to do some TSP based visualisations. In this guide, I’ll go through what I tried and the end results, in the hope that it may be a helpful reference for the future.

3.1. Generating the problem

using Random, Metaheuristics, WGLMakie, Bonito
Page(exportable=true, offline=true)

Here is a quick way to generate some data for which one can try to find the best (i.e. lowest cost/shortest) Hamiltonian cycle. This was taken from this JuMP tutorial.

function generate_distance_matrix(n; random_seed = 1)
    rng = Random.MersenneTwister(random_seed)
    X_coord = 100 * rand(rng, n)
    Y_coord = 100 * rand(rng, n)
    d = [sqrt((X_coord[i] - X_coord[j])^2 + (Y_coord[i] - Y_coord[j])^2) for i in 1:n, j in 1:n]
    return d, X_coord, Y_coord
end

Let’s pick a specific \(n\) and observe what the graph looks like.

n = 40
d, X_coord, Y_coord = generate_distance_matrix(n)

scatter(X_coord, Y_coord)

3.2. Solving the problem

When it comes to evolutionary algorithms, I have encountered 5 Julia packages: Metaheuristics.jl, Evolutionary.jl, Optim.jl, GeneticAlgorithms.jl, and GAFramework.jl. This is not meant to be an exhaustive list and thus may be missing other packages.

Of these 5, the last two are specifically focused on GAs and the third doesn’t have GAs but other algorithms. I’ve decided to try out the first two only, in part because they are both wrapped by Optimization.jl, which I thought may make plugging one out and the other in simpler.

3.2.1. Evolutionary.jl

Having decided to try Evolutionary.jl first for no reason, I had to formulate and feed the problem in a form acceptable by the package. Here, I ran into a problem: its documentation states that Evolutionary accepts constraints of the form

\[ \begin{align}\begin{aligned}l_x\leq x\leq u_x\\l_c \leq c(x)\leq u_c\end{aligned}\end{align} \]

(though presumably, those should be \(x_j\) instead of \(x\).) This on its own is not a problem, for example one can use the Dantzig-Fulkerson-Johnson formulation of TSP as an integer linear program, which would satisfy this form. But I couldn’t figure out the Evolutionary interface to supply the function \(c\).

Since it has a wrapper, I tried using Optimization.jl as an interface, but I still kept running into issues, and thus gave up.

3.2.2. Metaheuristics.jl

This ended up being much easier, primarily due to baked in support for SearchSpaces.jl (by the same author) and the provided PermutationSpace search space, offering a much more intuitive formulation of the TSP problem.

bounds = PermutationSpace(n)

The permutation representation gives the order of visiting various locations, and automatically incorporates constraints such as visiting every city and only once and without subloops. This way, the objective function also becomes very easy to define.

function tsp_distance(x)
    total = 0
    for i in 1:(n-1)
        total += d[x[i],x[i+1]]
    end
    total += d[x[n],x[1]]
    return total
end

One can then specify settings like iteration limit, tolerances and more.

options = Options(iterations=1000, seed=1)

Options
=======
  rng:             TaskLocalRNG()
  seed:            1
  x_tol:           1.0e-8
  f_tol:           1.0e-12
  g_tol:           0.0
  h_tol:           0.0
  debug:           false
  verbose:         false
  f_tol_rel:       2.220446049250313e-16
  time_limit:      Inf
  iterations:      1000
  f_calls_limit:   0.0
  store_convergence: false
  parallel_evaluation: false

Then we need to pick the algorithm to use. I’ll use the genetic algorithm implementation but the package offers others too, such as particle swarm.

ga = GA(;p_mutation = 0.5,
        p_crossover = 0.9,
        initializer = RandomPermutation(N=100), 
        crossover=OrderCrossover(), 
        mutation=SlightMutation(),
        options=options);

The last step is to run the optimization.

result = optimize(tsp_distance, bounds, ga)

Optimization Result
===================
  Iteration:       587
  Minimum:         579.164
  Minimizer:       [24, 1, 38, …, 9]
  Function calls:  58700
  Total time:      0.7615 s
  Stop reason:     Due to Convergence Termination criterion.

With this at hand, we can access the best result obtained and see what it looks like.

sol = minimizer(result)
fig, ax, splot = scatter(X_coord, Y_coord)
lines!(ax, X_coord[sol], Y_coord[sol])  # reorder vector using permutation

endpoints = sol[[1,n]]
lines!(ax, X_coord[endpoints], Y_coord[endpoints])  # connect the cycle
fig

3.3. Animations

Animations in Makie work by creating plots where every changing element is an Observable, and changing the values of these Observables to obtain different frames.

Suppose we want to animate the above plot to show the solutions obtained in each iteration. One fundamental element that is changing across iterations is the best solution for a given generation. We will keep track of it with an Observable.

perm = Observable(collect(1:n))

Observable([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  31, 32, 33, 34, 35, 36, 37, 38, 39, 40])

Here, collect(1:n) acts as a starting value. The type of Observables are fixed when they are created, so it can be important to specify that, but in our case, every permutation is represented as Vector{Int32}, so we are fine.

We can also keep track of the iteration to display in each frame.

it = Observable(0)

Observable(0)

Similar to the previous plots on this guide, we first plot the scatter plot.

fig, ax, splot = scatter(X_coord, Y_coord, axis=(title= @lift("t=$($it)"), ))

Here, the @lift macro creates a new Observable that depends on an already existing one. This way, when we change the value of it for example to \(1\), Makie will now that the title needs to become t=1. You can read more about this here, but the important part is that they are now linked, and that the syntax requires the Observable to be prepended with a dollar sign (which is inside a $(..) since we are also doing string formatting).

To keep things a bit tidier, I define a couple more things as linked Observables.

# Lines reordered according to the new permutation
X_lines = @lift(X_coord[to_value($perm)])
Y_lines = @lift(Y_coord[to_value($perm)])

# Endpoints for the new permutation
X_endpoints = @lift(X_coord[to_value($perm)[[1,n]]])
Y_endpoints = @lift(Y_coord[to_value($perm)[[1,n]]])

In the above, since perm is an Observable, one cannot index a Vector by it. to_value is used to access the permutation perm contains, which results in a valid reordering operation. Using these newly defined Observables, we can plot the starting point of the graph.

lines!(ax, X_lines, Y_lines)
lines!(ax, X_endpoints, Y_endpoints)
fig

Now, the basic examples of making animations in Makie follow the form

record(fig, "filename.gif", iterator) do it
    ...
end

or (an equivalent way of writing the same signature)

record(func, fig, "filename.gif", iterator)

where func is called for every iteration of iterator to create a new frame to be recorded.

It is not clear to me how to use this signature, since we are not going iteration by iteration; after calling optimize, Metaheuristics does everything for us. However, reading the documentation for record reveals an alternative:

record(func, fig, path)

where the iteration and frame creation responsibilities are manually handled. We can combine this with the logger mechanism of optimize, which allows one to provide a function that will be called with the current State after each iteration.

function logger(io, state)
    # get current best solution
    sol = minimizer(state)
    # update the Observables
    perm[] = sol
    it[] = state.iteration
    # record new frame of animation
    recordframe!(io)
end

Ultimately, we can obtain the animation with

record(fig, "tsp.mp4") do io
    optimize(tsp_distance, bounds, ga, 
        logger=(st)->logger(io, st)  # optimize only provides the state, so we need to pass everything properly
    )
end

3.4. Interactive Example

We can use Bonito.jl to obtain interactive examples. Since this website is static and does not communicate with an underlying Julia process, all the interactivity needs to happen in HTML/JavaScript. The record_states function records the possible states of a Widget and how it affects related Observables, as an alternative to programming the behavior in JavaScript.

Warning

The docs warn that record_states is experimental and unoptimized. So you may want to be careful if you’ll use it a lot.

To give a specific example, the HTML files for other pages are about ~30kB, but this one is about ~30MB.

First, the best permutations are reobtained, because I didn’t save them previously.

best_perms = []
function logger(st)
    push!(best_perms, minimizer(st))
end

ga = GA(;p_mutation = 0.5,
        p_crossover = 0.9,
        initializer = RandomPermutation(N=100), 
        crossover=OrderCrossover(), 
        mutation=SlightMutation(),
        options=options);

result = optimize(tsp_distance, bounds, ga, logger=logger)

Optimization Result
===================
  Iteration:       587
  Minimum:         579.164
  Minimizer:       [24, 1, 38, …, 9]
  Function calls:  58700
  Total time:      0.8494 s
  Stop reason:     Due to Convergence Termination criterion.

App() do session::Session
    it_slider = Slider(1:result.iteration)
    lab = @lift(string("t=", $(it_slider.value)))
    perm = map((x) -> best_perms[x], it_slider)
    fig, ax, splot = scatter(X_coord, Y_coord, axis=(title=lab, ))
    X_lines = @lift(X_coord[to_value($perm)])
    Y_lines = @lift(Y_coord[to_value($perm)])
    X_endpoints = @lift(X_coord[to_value($perm)[[1,n]]])
    Y_endpoints = @lift(Y_coord[to_value($perm)[[1,n]]])
    lines!(X_lines, Y_lines)
    lines!(X_endpoints, Y_endpoints)
    return Bonito.record_states(session, DOM.div(it_slider, it_slider.value, fig))
end

1