Linear Programming With A Childs Puzzle

Around 15 years ago I received this puzzle for Christmas, it had a simple goal; put all the tiles in place such that the bug's tops and bottoms and match. There are only 9 tiles and this was a children's toy, surely I'd have it done in a day or two...

But I never did and if it hadn't come solved in the box, I'd have never believed it was possible. Every few years I'd find it and have a go for a few hours then give up.

I'd learn something new at school and later university that would help but it would never quite be enough. When I came across this puzzle most recently two things were different; I had ample time due to lockdown, and I'd been learning python for the last few years.

Encoding The Puzzle

As you can see above, as well as the 9 tiles, there are 4 types of bug, each with a top and bottom, so the first step was to encode this puzzle in such a way that I could try combinations of tiles and determine if this was a valid solution or not.

I did this by giving each bug a numerical value, where the top was positive and the bottom negative, that way if the pairs were from the same bug, they would sum to 0.

Cricket Head = +1
Cricket Tail = -1
Grub Head = +2
Grub Tail = -2
Ladybird Head = +3
Ladybird Tail = -3
Spider Head = +4

Spider Tail = -4

As a tile each contains 4 bug sections but the ordering of these sections is important, I decided to encode each tile as a 4 value numpy array where the first value is the right most segment and goes around the tile clockwise.

So the tile above would be expressed using:

piece_8 = np.array([3, -3, 1, 4])

With this, we can take each vector as a row and form a matrix of size 9x4 using the stack function in numpy. This matrix would correspond in the puzzle to placing the 9 tiles in some combination of positions, each tile being oriented in one of four ways.

The next step from here was to create a function that could determine if any combination of pieces was a valid solution or not.

As you can see, there are 12 bugs that must match and these are always in the same places, so we can determine that for a given matrix, there are 24 specific elements that form 12 pairs, where each pair must sum to 0. Thus we can now determine if a matrix is a solution or not.

The Brute Force Approach

My first idea was to simply try every combination of tile placement and orientation until a solution was found, but for this we would need functions for placing and rotating each tile.

As each of the 9 tiles can be placed in 4 ways due to its orientation, independent of all the others, we find that one setting of the 9 tiles would result in 4^9 = 262,144 matrices to check for a valid solution.

On top of this, there are 9! = 362,880 different ways to place the tiles, this would lead to a maximum of 95,126,814,720 matrices to check.

Upon timing my checker function, I found it could check a matrix in a fraction of a second, but even if this was brought down to say 10 milliseconds with efficiency improvements, this would still lead to a maximal search time of approximately 30 years!

A more elegant solution was needed.

A Linear Programming Approach

When discussing this problem with a friend, I was pointed in the direction of Linear Programming, an area of mathematical optimization that is concerned with finding the optimal solution to a problem where the constraints have linear relationships to each other.

My favourite example of this being a Sudoku puzzle, where using a brute force approach would involve a huge number of calculations but with linear programming, can be solved almost instantly.

The more I read in to the topic, the more I realised that this was exactly the sort of problem I had. A quick search found a python library named "constraint" that would allow me to encode the vectors I had as a constraint problem. This led to the following function:

import constraint
import numpy as np
import time




def constraint_checker(piece_1, piece_2, piece_3, piece_4, piece_5,
                       piece_6, piece_7, piece_8, piece_9, timing=False):


    # Find the start time if we want the execution time for this function
    if timing:
        start_time = time.time()


    # Create the set of rotations for each piece
    piece_1_set = [list(np.roll(piece_1, i)) for i in range(4)]
    piece_2_set = [list(np.roll(piece_2, i)) for i in range(4)]
    piece_3_set = [list(np.roll(piece_3, i)) for i in range(4)]
    piece_4_set = [list(np.roll(piece_4, i)) for i in range(4)]
    piece_5_set = [list(np.roll(piece_5, i)) for i in range(4)]
    piece_6_set = [list(np.roll(piece_6, i)) for i in range(4)]
    piece_7_set = [list(np.roll(piece_7, i)) for i in range(4)]
    piece_8_set = [list(np.roll(piece_8, i)) for i in range(4)]
    piece_9_set = [list(np.roll(piece_9, i)) for i in range(4)]


    # Initialise the constraint problem class
    problem = constraint.Problem()


    # Add the variables to the problem and their domain sets
    problem.addVariable("1", piece_1_set)
    problem.addVariable("2", piece_2_set)
    problem.addVariable("3", piece_3_set)
    problem.addVariable("4", piece_4_set)
    problem.addVariable("5", piece_5_set)
    problem.addVariable("6", piece_6_set)
    problem.addVariable("7", piece_7_set)
    problem.addVariable("8", piece_8_set)
    problem.addVariable("9", piece_9_set)


    # Add the constraints required to solve the puzzle
    problem.addConstraint(lambda i, j: i[1] + j[3] == 0, ("1", "2"))
    problem.addConstraint(lambda i, j: i[1] + j[3] == 0, ("4", "5"))
    problem.addConstraint(lambda i, j: i[1] + j[3] == 0, ("7", "8"))
    problem.addConstraint(lambda i, j: i[1] + j[3] == 0, ("2", "3"))
    problem.addConstraint(lambda i, j: i[1] + j[3] == 0, ("5", "6"))
    problem.addConstraint(lambda i, j: i[1] + j[3] == 0, ("8", "9"))
    problem.addConstraint(lambda i, j: i[2] + j[0] == 0, ("1", "4"))
    problem.addConstraint(lambda i, j: i[2] + j[0] == 0, ("4", "7"))
    problem.addConstraint(lambda i, j: i[2] + j[0] == 0, ("2", "5"))
    problem.addConstraint(lambda i, j: i[2] + j[0] == 0, ("5", "8"))
    problem.addConstraint(lambda i, j: i[2] + j[0] == 0, ("3", "6"))
    problem.addConstraint(lambda i, j: i[2] + j[0] == 0, ("6", "9"))


    # Give a list of all solutions
    solutions = problem.getSolutions()


    # Find end time and difference, return this with solutions
    if timing:
        end_time = time.time()
        execution_time = end_time - start_time
        return solutions, execution_time
    else:

        
        return solutions

The above function would allow me to check one of the 9! combinations of tile placements in less than 10 milliseconds, a much greater improvement and small enough that when looped over each permutation, a solution could be found within an hour.

I wrote the rest of the script and set my laptop to work, 10 minutes later I had a possible solution! Once I had translated the vector encoding back to the physical pieces I was pleased to see I had in fact brought this matter to a close as shown below.

One thing that I did wonder about was how I found the solution in less than one sixth of the maximum time, but when placing the tiles on the board, I noticed the whole board could be rotated 4 times, so in effect, in my regime there would in fact be 4 unique solutions to this problem for each valid "human" solution. This would make my 10 minute time seem more reasonable within a 15 minute maximum.

In conclusion, although all of this was really just a bit of fun to kill some time during lockdown, to me it shows how even in the simplest of places, we can find unexpected layers of complexity and never even notice them until we come back armed new techniques and skills.

All code available at on my GitHub and any questions or comments on this or any other interesting/obscure puzzles, feel free to connect over LinkedIn, I'm always looking to learn something and bring the interesting to the mundane.