Julia Community 🟣: Evan Wright

MILP with JuMP and HiGHS

Evan Wright — Fri, 17 Feb 2023 05:33:52 +0000

In my previous post, I mentioned that the problem (Advent of Code 2018 day 23) can be reformulated as a mixed-integer linear program (MILP). In this post, we'll walk through a solution using JuMP.jl and HiGHS.jl. The formulation is based on this Reddit comment.

Input parsing is the same as last time. We set up the JuMP problem by defining variables x, y, and z as integers and the vector variable botsinrange as binary. This time, the value of botsinrange[i] is 1 if the point (x, y, z) is in range of bot i, otherwise 0.

using JuMP
using HiGHS

struct Nanobot
    x::Int
    y::Int
    z::Int
    radius::Int
end
function processinput(filename)
    puzin = [[parse(Int, m.match) for m in eachmatch(r"(\-?\d+)", line)]
             for line in readlines(filename)]
    bots = [Nanobot(l...) for l in puzin]
    return bots
end
bots = processinput("input201823.txt")

model = Model(HiGHS.Optimizer)
set_silent(model)
@variable(model, x, Int)
@variable(model, y, Int)
@variable(model, z, Int)
@variable(model, botsinrange[1:length(bots)], Bin)

Next, define constraints of the form

abs(b.x - x) + abs(b.y - y) + abs(b.z - z) <= b.radius + (1 - botsinrange[i])*slack

where i is the index of bot b. For the moment, assume slack is infinity. Then if we think bot i is out of range (botsinrange[i] = 0), the constraint holds (right hand side is Inf). Otherwise, if we think bot i is in range (botsinrange[i] = 1), the constraint holds if and only if abs(b.x - x) + abs(b.y - y) + abs(b.z - z) <= b.radius, i.e. the bot has to actually be in range. We will maximize the sum of botsinrange, so the optimizer will try to set botsinrange[i] = 1 for as many i as possible, and in those cases, it must also ensure the corresponding constraint holds.

So far so good, but there are two minor complications. First, the abs function is non-linear, so we can't have it in the constraint definition. However, note that if -x < c and x < c, then abs(x) < c. If x + y < c and -x + y < c and x - y < c and -x - y < c, then abs(x) + abs(y) < c. A similar argument holds for 3 variables, so we just need to create several linear constraints that imply the constraint defined above.

Second, we can't actually use Inf for the slack for numerical stability reasons; we'll have to use a large finite number. The slack can't be too low because we need the constraint to be relaxed if botsinrange[i] = 0. Consider the worst case, where the left hand side of the constraint is as large as possible. This may occur when (x, y, z) is at the "opposite" corner of 3d space from a bot i. The largest coordinate or radius in our puzzle input is

maxcoord = maximum(b -> max(abs(b.x), abs(b.y), abs(b.z), b.radius), bots)

In the worst case, we are comparing a point like (maxcoord, maxcoord, maxcoord) to (-maxcoord, -maxcoord, -maxcoord), so the distance between those points is 2*3*maxcoord. We have additional distance for the bots' sensor radius, so we may use slack = 2*4*maxcoord.

Putting it all together, our constraints are

slack = 2*4*maximum(b -> max(abs(b.x), abs(b.y), abs(b.z), b.radius), bots)
for (i,b) in enumerate(bots)
    for absmult in Iterators.product((-1,1),(-1,1),(-1,1))
        @constraint(model, absmult[1]*(b.x - x) + absmult[2]*(b.y - y) +
                            absmult[3]*(b.z - z) <= b.radius + (1-botsinrange[i])*slack)
    end
end

Finally, we can optimize the model.

@objective(model, Max, sum(botsinrange))
optimize!(model)
println("""x=$(round(Int, value(x))), y=$(round(Int, value(y))), z=$(round(Int, value(z)))
ans=$(round(Int, value(x)) + round(Int, value(y)) + round(Int, value(z)))""")

For my input, we find the correct answer, but recall the problem statement asked if there are multiple points with maximum nanobot coverage, then we should choose the point closest to the origin. So, we can add a constraint that the objective function must be at least as high as we just found¹

function coordisinrange(x, y, z, bot::Nanobot)
    return abs(bot.x - x) + abs(bot.y - y) + abs(bot.z - z) <= bot.radius
end
maxbots = count(b -> coordisinrange(round(Int, value(x)),
         round(Int, value(y)), round(Int, value(z)), b), bots)
@constraint(model, maxbots <= sum(botsinrange))

and add a penalty term to our objective function to minimize abs(x) + abs(y) + abs(z).² We have the same non-linear problem with abs as before, so we have to use a similar trick. The objective function becomes sum(botsinrange) - 0.00000001 * (xt + yt + zt) where the penalty is set by trial and error until the optimizer finds a solution. The coordinates are 8 digit numbers, so the penalty should be similar in scale to 1 unit of the objective---i.e. having one more bot in range.

@variable(model, xt, Int)
@variable(model, yt, Int)
@variable(model, zt, Int)
@constraint(model, x <= xt); @constraint(model, -x <= xt);
@constraint(model, y <= yt); @constraint(model, -y <= yt);
@constraint(model, z <= zt); @constraint(model, -z <= zt);

@constraint(model, maxbots <= sum(botsinrange))
@objective(model, Max, sum(botsinrange) - 0.00000001 * (xt + yt + zt))
optimize!(model)
println("""x=$(round(Int, value(x))), y=$(round(Int, value(y))), z=$(round(Int, value(z)))
ans=$(round(Int, value(x)) + round(Int, value(y)) + round(Int, value(z)))""")

The solution is the same answer as before---which we know is correct thanks to Z3---but I'm curious if any other input needs the tie-breaking rule.

Discussion

The MILP solution feels a bit hacky compared to Z3, but it is faster at around 0.3-1.5 seconds without compilation, compared to 50+ seconds for Z3. It took trial and error to get the right slack, and even if the slack is high enough, the solver may provide an incorrect solution.³ For my input, slack of 1e11 results in an incorrect solution, while 1e10, 1e12-1e14 are ok.

Handling the L1 norm (abs) is annoying as well. I tried to use MathOptInterface.NormOneCone,

for (i,b) in enumerate(bots)
    @constraint(model2,
     [b.radius + (1-botsinrange[i])*Int(1e10); x - [b.x, b.y, b.z]] in MOI.NormOneCone(4))
end

but there is a performance issue. With constraints from only 280 bots (the problem has 1000), the solver takes about 4 seconds, but with 290 bots it hangs (or takes longer than I care to wait).

In summary, MILP solvers can solve this problem more quickly than Z3 with a minor reformulation. There may be a better formulation that is less finicky. If you know of one, let me know.

Full Julia script

using JuMP
using HiGHS

struct Nanobot
    x::Int
    y::Int
    z::Int
    radius::Int
end

function processinput(filename)
    puzin = [[parse(Int, m.match) for m in eachmatch(r"(\-?\d+)", line)]
             for line in readlines(filename)]
    bots = [Nanobot(l...) for l in puzin]
    return bots
end

function coordisinrange(x, y, z, bot::Nanobot)
    return abs(bot.x - x) + abs(bot.y - y) + abs(bot.z - z) <= bot.radius
end

function main(filename)
    bots = processinput(filename)
    model = Model(HiGHS.Optimizer)
    slack = 2*4*maximum(b -> max(abs(b.x), abs(b.y), abs(b.z), b.radius), bots)
    set_silent(model)
    @variable(model, x, Int)
    @variable(model, y, Int)
    @variable(model, z, Int)
    @variable(model, botsinrange[1:length(bots)], Bin)
    for (i,b) in enumerate(bots)
        for absmult in Iterators.product((-1,1),(-1,1),(-1,1))
            @constraint(model, absmult[1]*(b.x - x) + absmult[2]*(b.y - y) +
                                absmult[3]*(b.z - z) <= b.radius + (1-botsinrange[i])*slack)
        end
    end
    @objective(model, Max, sum(botsinrange))
    optimize!(model)
    println("Step 1 solution:")
    println("""x=$(round(Int, value(x))), y=$(round(Int, value(y))), z=$(round(Int, value(z))),
     ans=$(round(Int, value(x)) + round(Int, value(y)) + round(Int, value(z)))""")

    # Step 2, with lower bound constraint for number of bots in range.
    maxbots = count(b -> coordisinrange(round(Int, value(x)),
         round(Int, value(y)), round(Int, value(z)), b), bots)
    @constraint(model, maxbots <= sum(botsinrange))

    @variable(model, xt, Int)
    @variable(model, yt, Int)
    @variable(model, zt, Int)
    @constraint(model, x <= xt); @constraint(model, -x <= xt);
    @constraint(model, y <= yt); @constraint(model, -y <= yt);
    @constraint(model, z <= zt); @constraint(model, -z <= zt);

    @objective(model, Max, sum(botsinrange) - 0.00000001 * (xt + yt + zt))
    optimize!(model)
    println("Step 2 solution:")
    println("""x=$(round(Int, value(x))), y=$(round(Int, value(y))), z=$(round(Int, value(z)))
     ans=$(round(Int, value(x)) + round(Int, value(y)) + round(Int, value(z)))""")

    return nothing
end

@time main("input201823.txt")

Output (second run to avoid compilation time):

Step 1 solution:
x=44166763, y=43480550, z=38585775,
ans=126233088
Step 2 solution:
x=44166763, y=43480550, z=38585775
ans=126233088
  1.688667 seconds (590.15 k allocations: 32.233 MiB, 0.86% gc time)

You may be wondering why I didn't use round(Int, objective_value(model)). I found that for some values of slack, the objective_value is not the value of the objective at the solution (x, y, z). For example, the solution for step one results in an objective_value of 969 when the true number of bots in range of (x, y, z) is 972. 969 is still larger than the objective function evaluated at any other point, though. There may be a setting for HiGHS---that I'm unaware of---to resolve this issue. ↩
We could do this in the first step as well, but I'm more confident in finding the global maxbots without the penalty in the first step, then trying to minimize the L1 norm of the solution. ↩
If the slack is too low, then the problem is incorrectly specified. In the extreme case where slack is 0, the problem is to find a point within every bot's range, which isn't feasible. However, there shouldn't be an issue---mathematically---with a slack that is too large. ↩

Z3.jl - Optimization Examples Solved with Julia

Evan Wright — Wed, 15 Feb 2023 05:37:28 +0000

Z3 is a satisfiability modulo theories (SMT) solver callable from Julia with Z3.jl. Documentation is somewhat lacking, so I'd like to share a worked example of optimizing an integer programming problem with Z3.jl.

We will solve Advent of Code 2018 day 23, so if you'd like to avoid spoilers, stop reading. The problem provides a list of nanobots defined by a center point in 3d space and radius, measured by Manhattan distance. For part 2 of the problem, we must find the (integer) point that lies within range of the largest number of nanobots. If there are multiple such points, we must find the one closest to the origin (0, 0, 0).

First, we'll define a struct to hold the nanobot data and process the input
file.

using Z3

struct Nanobot
    x::Int
    y::Int
    z::Int
    radius::Int
end

function processinput(filename)
    puzin = [[parse(Int, m.match) for m in eachmatch(r"(\-?\d+)", line)]
             for line in readlines(filename)]
    bots = [Nanobot(l...) for l in puzin]
    return bots
end

bots = processinput("input201823.txt")

You can use the following example input if you don't want to log in to Advent of Code.

pos=<10,12,12>, r=2
pos=<12,14,12>, r=2
pos=<16,12,12>, r=4
pos=<14,14,14>, r=6
pos=<50,50,50>, r=200
pos=<10,10,10>, r=5

Let's think about how we may define the problem as an integer optimization that Z3 can understand. Consider a potential solution point (x, y, z). The point is outside the radius of nanobot bot if and only if abs(bot.x - x) + abs(bot.y - y) + abs(bot.z - z) > bot.radius. We want to minimize the number bots for which that statement is true. We can collect the statements in an array, and sum the boolean values (1 = true). The resulting sum (call it noutofrange) will be our objective to minimize.

Finally, conditional on minimizing noutofrange, we want to choose the point closest to the origin, or minimize abs(x) + abs(y) + abs(z) (call this value cost). Z3 uses by default a lexicographic priority of objectives. It solves first for the objective that is declared first, so this should be easy.

Define the Z3 context and necessary variables. We'll see why we need z3zero and z3one in a moment.

ctx = Context()
x = int_const(ctx, "x")
y = int_const(ctx, "y")
z = int_const(ctx, "z")
z3zero = int_val(ctx, 0)
z3one = int_val(ctx, 1)
noutofrange = int_const(ctx, "noutofrange")
cost = int_const(ctx, "cost")

The array of boolean values (converted to integer 1 or 0) measuring whether the point is within range of each bot may now be defined as

botsoutofrange = [ite(abs(b.x - x) + abs(b.y - y) + abs(b.z - z) > b.radius,
                      z3one, z3zero) for b in bots]

ite(a, b, c) (if then else) is a Z3 method that means

if a
    return b
else
    return c
end

We need b and c to be Z3 expressions, hence the use of our previously defined z3one and z3zero.

Now, we initialize the optimization problem with the "constraints" that define our objective variables,

opt = Optimize(ctx)
add(opt, noutofrange == sum(botsoutofrange))
add(opt, cost == abs(x) + abs(y) + abs(z))

and finally optimize each of the objectives in order.

minimize(opt, noutofrange)
minimize(opt, cost)
res = check(opt)

Confirm that the problem is satisfiable and print the optimal point.

@assert res == Z3.sat
m = get_model(opt)
for (k, v) in consts(m)
    println("$k = $v")
end

For my input, this prints

noutofrange = 28
cost = 126233088
z = 38585775
y = 43480550
x = 44166763

It takes about a minute to solve on my PC. That may seem slow, but the human time required to code a faster algorithm is a lot longer than a minute.¹ Z3 is one of those tools that isn't often applicable, but when it is, it feels like magic.

Full Julia script

using Z3

struct Nanobot
    x::Int
    y::Int
    z::Int
    radius::Int
end

function processinput(filename)
    puzin = [[parse(Int, m.match) for m in eachmatch(r"(\-?\d+)", line)]
             for line in readlines(filename)]
    bots = [Nanobot(l...) for l in puzin]
    return bots
end

function main(filename)
    bots = processinput(filename)
    ctx = Context()
    x = int_const(ctx, "x")
    y = int_const(ctx, "y")
    z = int_const(ctx, "z")
    z3zero = int_val(ctx, 0)
    z3one = int_val(ctx, 1)
    noutofrange = int_const(ctx, "noutofrange")
    cost = int_const(ctx, "cost")
    botsoutofrange = [ite(abs(b.x - x) + abs(b.y - y) + abs(b.z - z) > b.radius, z3one, z3zero)
                      for b in bots]
    opt = Optimize(ctx)
    add(opt, noutofrange == sum(botsoutofrange))
    add(opt, cost == abs(x) + abs(y) + abs(z))
    minimize(opt, noutofrange)
    minimize(opt, cost)
    res = check(opt)
    @assert res == Z3.sat
    m = get_model(opt)
    for (k, v) in consts(m)
        println("$k = $v")
    end
    return nothing
end

@time main("input201823.txt")

Output:

noutofrange = 28
cost = 126233088
z = 38585775
y = 43480550
x = 44166763
 52.030344 seconds (71.38 k allocations: 3.053 MiB, 0.07% compilation time)

Comparison with Python

The script below is the same problem solved with the Python package z3-solver. I had to manually define abs to work with Z3 in Python, but Python didn't require ite or definition of the 0 and 1 Z3 expressions like Julia. I also had occasional segfaults in the Julia version, although it wasn't reproducible when executing the file in a fresh REPL. It may be related to this issue: https://github.com/ahumenberger/Z3.jl/issues/12.

import z3
import re
import time

with open('input/input201823.txt', 'r') as f:
    actualinput = f.read()

bots = [list(map(int, re.findall('(\-?\d+)', row))) for row in actualinput.splitlines()]

# https://stackoverflow.com/questions/22547988/how-to-calculate-absolute-value-in-z3-or-z3py
def z3abs(x):
    return z3.If(x >= 0,x,-x)

x = z3.Int('x')
y = z3.Int('y')
z = z3.Int('z')
noutofrange = z3.Int('noutofrange')
cost = z3.Int('cost')

botsoutofrange = [z3abs(b[0]-x) + z3abs(b[1]-y) + z3abs(b[2]-z) > b[3] for b in bots]

opt = z3.Optimize()
opt.add(noutofrange == z3.Sum(botsoutofrange))
opt.add(cost == z3abs(x) + z3abs(y) + z3abs(z))
# Z3 uses by default a lexicographic priority of objectives. It solves first for the objective that is declared first.
lowestoutofrange = opt.minimize(noutofrange)
closesttoorigin = opt.minimize(cost)
print("Checking optimization...")
start_time = time.time()
opt.check()
print(f"Found solution in {time.time() - start_time:.2f}s")
print(opt.model())

The problem can be reformulated as a mixed-integer linear program, which is solvable nearly instantly by an optimizer like HiGHS, but that still requires some human effort to think through and implement the reformulation. ↩

Solving Peg Solitaire with Julia

Evan Wright — Sun, 19 Jun 2022 03:51:07 +0000

Peg solitaire is a singleplayer board game with the objective to remove all game pieces (pegs or marbles) except one from the board by "jumping" them with another peg. The 15-hole triangular variant is commonly found in Cracker Barrel restaurants in the US. We will solve this variant with Julia.

The triangular board is laid out as a hexagonal grid, so jumps can occur by moving northwest, northeast, east, southeast, southwest, or west. To simplify our code, we will represent the board as a BitMatrix, a 2d array of 1s and 0s, and restrict moves to valid indices.

INDICES = Vector{Tuple{Int, Int}}()
board = falses(5, 9)
for col in axes(board, 2)
    for row in axes(board, 1)
        if ((row % 2 == col % 2) && (col >= row) &&
                (col <= size(board, 2) + 1 - row))
            board[row, col] = true
            push!(INDICES, (row, col))
        end
    end
end
println("The board is represented as")
display(board)
println("Valid indices are")
println(INDICES)

The board is represented as
5×9 BitMatrix:
 1  0  1  0  1  0  1  0  1
 0  1  0  1  0  1  0  1  0
 0  0  1  0  1  0  1  0  0
 0  0  0  1  0  1  0  0  0
 0  0  0  0  1  0  0  0  0
Valid indices are
[(1, 1), (2, 2), (1, 3), (3, 3), (2, 4), (4, 4),
(1, 5), (3, 5), (5, 5), (2, 6), (4, 6), (1, 7),
(3, 7), (2, 8), (1, 9)]

To make our life a bit easier, let's first write a function to render the board. I love printstyled by the way. If you use the Julia REPL or a supported terminal, the pegs will be cyan.

function showboard(board; upsidedown=false)
    if upsidedown
        iiter = reverse(axes(board, 1))
    else
        iiter = axes(board, 1)
    end
    println("\n------------------")
    for i in iiter
        for j in axes(board, 2)
            if board[i, j]
                printstyled('█'; color=:cyan)
            elseif (i % 2 == j % 2) && (j >= i) &&
                    (j <= size(board, 2) + 1 - i)
                printstyled('o'; color=:nothing)
            else
                printstyled('░'; color=:nothing)
            end
        end
        if i == first(iiter) + (last(iiter) - first(iiter))÷2
            print(" pegs: $(sum(board))")
        end
        print("\n")
    end
    print("------------------\n")
end

board[1, 1] = false # remove a peg for illustration
showboard(board)

------------------
o░█░█░█░█
░█░█░█░█░
░░█░█░█░░ pegs: 14
░░░█░█░░░
░░░░█░░░░
------------------

This board is a potential starting position. All holes are filled with pegs except the top left one.

Next, we need to determine the potential moves. In other words, given the current state of the board, what are the potential states of the board after one jump? For each peg, consider its neighbor positions (northwest, northeast, east, southeast, southwest, or west). The neighbor position must contain a peg, so we have something to jump. In addition, the neighbor's neighbor in the same direction must be empty, so we have a hole to jump to.

Using the following figure, suppose i is the location of the peg that will jump. Then, locations marked n are its neighbors, and n2 are the neighbors' neighbors at a valid index. There are two potential moves: (1) jump from 3,5 to 1,3 and remove the peg at 2,4, and (2) jump from 3,5 to 1,7 and remove the peg at 2,6. As long as the applicable neighbor positions have a peg and the neighbors' neighbors are empty, these moves are valid.

  1  2  3  4  5  6  7  8  9 
1       n2          n2   
2          n     n    
3       n     i     n   
4          n     n    
5

We must repeat this process, checking for valid moves for each peg in the board. The following function returns a vector of BitMatrix representing board states reachable in one move from the current board state (an input BitMatrix). The input indices is a vector of tuples representing valid indices for peg locations, so we don't have to compute them for each call of the function. mdist is a helper function to find valid neighbor locations.

function mdist(a, b)
    return abs(a[1] - b[1]) + abs(a[2] - b[2])
end

function getnextstates(board, indices)
    nextstates = Vector{BitMatrix}()
    for ind in indices
        # (1) there must be a peg to use to jump
        if !board[ind...]
            continue
        end
        neighbors = (n for n in indices if ind[2] != n[2] &&
                                           mdist(ind, n) <= 2)
        for n in neighbors
            # (2) the neighbor location must contain a peg
            if !board[n...]
                continue
            end
            n2 = (n[1] - (ind[1] - n[1]), n[2] - (ind[2] - n[2]))
            # (3) the neighbor's neighbor in the same direction
            # must be a valid location and empty
            if n2 in indices && !board[n2...]
                # If (1), (2), and (3) are satisfied,
                # make a copy of the current board.
                newboard = copy(board)
                # Remove the jumping peg from its current location.
                newboard[ind...] = false
                # Remove the jumped peg.
                newboard[n...] = false
                # Place the jumping peg in its new location.
                newboard[n2...] = true
                # Save the new board state.
                push!(nextstates, newboard)
            end
        end
    end
    return nextstates
end

Now that we have a way to compute the states reachable from the current state, we can apply depth-first search (DFS) to find a solution. A full explanation of DFS is outside the scope of this note. However, we use the iterative form (rather than recursive), and the general flow for this problem is as follows:

Initialize a vector of states to search (states), beginning with the current state (board). Initialize an empty vector of states representing the path to the current state (pathtocurrentstate).¹
Pop a state from the vector of states to search, and add it to the path.
Check if the state is a winner, i.e. the number of pegs is 1. If so, return the path.
Generate the next possible states, and append them to the vector to be searched next.
Go to 2.

function findsolution(board, indices)
    states = [copy(board)]
    pathtocurrentstate = [copy(board)]
    pop!(pathtocurrentstate)
    while !isempty(states)
        currentstate = pop!(states)
        while (length(pathtocurrentstate) > 0 &&
               sum(pathtocurrentstate[end]) <= sum(currentstate))
            # Use the sum as a measure of depth.
            # Remove any same or lower depth states before pushing
            # the current state to the path.
            # This is necessary when we recurse "up".
            pop!(pathtocurrentstate)
        end
        push!(pathtocurrentstate, currentstate)
        if sum(currentstate) <= 1
            # This is the winning state
            return pathtocurrentstate
        end
        nextstates = getnextstates(currentstate, indices)
        if !isempty(nextstates)
            # If there are moves available from current_state,
            # push them to states (recursing down)
            append!(states, nextstates)
        end
    end
    return nothing
end

Now that we have a way to find a solution, let's write a single function to tie everything together.

function main(printsolution=true)
    # Initialize board
    INDICES = Vector{Tuple{Int, Int}}()
    board = falses(5, 9)
    for col in axes(board, 2)
        for row in axes(board, 1)
            if ((row % 2 == col % 2) && (col >= row) &&
                (col <= size(board, 2) + 1 - row))
                board[row, col] = true
                push!(INDICES, (row, col))
            end
        end
    end

    board[1,1] = false
    solution =  findsolution(board, INDICES)
    if printsolution
        if isnothing(solution)
            println("We didn't find any solution.")
        else
            println("Found the following solution:")
            for state in solution
                showboard(state)
            end
        end
    end
    return solution
end

main()

Found the following solution:

------------------
o░█░█░█░█
░█░█░█░█░
░░█░█░█░░ pegs: 14
░░░█░█░░░
░░░░█░░░░
------------------

------------------
█░o░o░█░█
░█░█░█░█░
░░█░█░█░░ pegs: 13
░░░█░█░░░
░░░░█░░░░
------------------

...

------------------
o░o░o░o░o
░o░o░o░o░
░░o░o░o░░ pegs: 2
░░░o░█░░░
░░░░█░░░░
------------------

------------------
o░o░o░o░o
░o░o░o░o░
░░o░o░█░░ pegs: 1
░░░o░o░░░
░░░░o░░░░
------------------

You'll have to run the code yourself if you want the full solution 😀.

Our method is reasonably fast,

@time main(false)
  0.003810 seconds (32.42 k allocations: 1.786 MiB)

but there is quite some room for improvement, especially if we want to find all solutions. getnextstates copies a lot. We could use an integer to directly represent the board state. I believe this is what BitMatrix does under the hood, but it probably has some overhead. We could cache the states we've already searched, including transformations of the board, which is symmetric. For more discussion on peg solitaire solution techniques or additional variants, refer to this page maintained by George Bell.

The purpose of pathtocurrentstate = [copy(board)] then pop!(pathtocurrentstate) is just to ensure pathtocurrentstate has the same type as states. ↩