Einstein Puzzle Solver

python
Author

Blaine Buxton

Published

August 20, 2023

Puzzle and idea

I love number puzzles and bought one inspired by Einstein on holiday recently. The idea is simple. There’s 16 tiles with numbers that can be placed in a 4x4 grid. The goal is for each row, column, and diagonal to sum up to 264. The tiles can be flipped so if tile has “66” imprinted on it, it can also represent “99”.

Here’s an example layout of tiles.

---------------------
| 18 | 89 | 98 | 61 |
-----+----+----+-----
| 68 | 91 | 88 | 16 |
-----+----+----+-----
| 81 | 19 | 66 | 99 |
-----+----+----+-----
| 96 | 69 | 11 | 86 |
---------------------

Things quickly escalated when I got the idea to write a solver that would use genetic programming concepts. The algorithm is simple. Start with a population of all scrambled states, then pick the best layouts out of the population. Exchange a few tiles in each of samples from the the best to create a new population along with a few new scrambled states. Repeat until we get a layout of tiles that satisfies the goal.

Puzzle Domain

First things first. A domain to model the puzzle was needed.

Source
import itertools
import random
from typing import Iterable, Sequence

# Type to represent the layout of tiles
State = list[int]

# Map the values of a number when flipped
FLIP_MAP: dict[int, int] = {
    1: 1,
    6: 9,
    8: 8,
    9: 6,
}

# All the tiles represented
ALL_TILES: State = [
    66,
    66,
    19,
    19,
    86,
    86,
    89,
    89,
    16,
    16,
    18,
    18,
    88,
    96,
    69,
    11,
]

# The goal that each row, column, and diagonal should sum up to
GOAL = 264


def flip(tile: int, flip_map: dict[int, int] = FLIP_MAP) -> int:
    """
    Take a tile and flip it. So, 89 will become 68.
    Calling it twice should return the original tile. 89->68->89
    """
    assert tile > 9 and tile < 100
    ones = tile % 10
    tens = tile // 10
    return flip_map[ones] * 10 + flip_map[tens]


def add_rows(tiles: State) -> list[int]:
    """
    Return the sum of all rows
    """
    return [sum(tiles[0:4]), sum(tiles[4:8]), sum(tiles[8:12]), sum(tiles[12:16])]


def add_columns(tiles: State) -> list[int]:
    """
    Return the sum of all columns
    """
    return [
        sum(each for i, each in enumerate(tiles) if i % 4 == 0),
        sum(each for i, each in enumerate(tiles) if i % 4 == 1),
        sum(each for i, each in enumerate(tiles) if i % 4 == 2),
        sum(each for i, each in enumerate(tiles) if i % 4 == 3),
    ]


def add_diagonals(tiles: State) -> list[int]:
    """
    Return the sum of the two diagonals
    """
    return [
        tiles[0] + tiles[5] + tiles[10] + tiles[15],
        tiles[12] + tiles[9] + tiles[6] + tiles[3],
    ]


def score_iter(tiles: State) -> Iterable[int]:
    """
    Return a score for each row, column, and diagonal.
    The score is the distance away from the goal.
    """
    return (
        abs(each - GOAL)
        for each in itertools.chain(
            add_rows(tiles), add_columns(tiles), add_diagonals(tiles)
        )
    )


def score(tiles: State) -> float:
    "Return a sum of how far away from goal for each row, column, and diagonal. This score is divided by the goal."
    return sum(each / float(GOAL) for each in score_iter(tiles))


def score_breakdown(tiles: State) -> list[float]:
    "Return a list of all the scores"
    return list(score_iter(tiles))


def chance(probability: float) -> bool:
    "Return True with probability between 1 and 0"
    return random.random() < probability


def exchange(
    tiles: State,
    probability_to_flip: float = 0.5,
    probability_to_exchange: float = 0.95,
) -> State:
    "Return a new state of tiles after exchanging two and possibly flipping either"
    choices = random.sample(range(len(tiles)), k=2)
    new_tiles = tiles.copy()
    for index in choices:
        if chance(probability_to_flip):
            new_tiles[index] = flip(new_tiles[index])
    if chance(probability_to_exchange):
        new_tiles[choices[0]], new_tiles[choices[1]] = (
            new_tiles[choices[1]],
            new_tiles[choices[0]],
        )
    return new_tiles


def scramble(tiles: State, probability_to_flip: float = 0.5) -> State:
    "Scramble all of the tiles"
    result = random.sample(tiles, k=len(tiles))
    return [flip(each) if chance(probability_to_flip) else each for each in result]

Domain Tests

To make sure the domain is correct, some tests were needed.

Source
import ipytest

ipytest.autoconfig()
ipytest.clean()

import math


def test_flip():
    def assert_flip(input: int, expected: int):
        assert flip(input) == expected
        assert flip(expected) == input

    assert_flip(66, 99)
    assert_flip(19, 61)
    assert_flip(86, 98)
    assert_flip(89, 68)
    assert_flip(16, 91)
    assert_flip(18, 81)
    assert_flip(88, 88)
    assert_flip(96, 96)
    assert_flip(69, 69)
    assert_flip(11, 11)


def test_add_rows():
    tiles = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
    actual = add_rows(tiles)
    assert actual == [10, 26, 42, 58]


def test_add_columns():
    tiles = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
    actual = add_columns(tiles)
    assert actual == [28, 32, 36, 40]


def test_add_diagonals():
    tiles = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
    actual = add_diagonals(tiles)
    assert actual == [1 + 6 + 11 + 16, 4 + 7 + 10 + 13]


def test_score():
    actual = score(ALL_TILES)
    assert math.isclose(actual, 2.708, abs_tol=0.001)


def test_solution():
    actual = score([91, 89, 68, 16, 18, 66, 81, 99, 86, 98, 19, 61, 69, 11, 96, 88])
    assert math.isclose(actual, 0.0, abs_tol=0.001)


def test_solution_another():
    actual = score([18, 86, 69, 91, 99, 61, 88, 16, 81, 19, 96, 68, 66, 98, 11, 89])
    assert math.isclose(actual, 0.0, abs_tol=0.001)


def test_score_breakdown():
    actual = score_breakdown(ALL_TILES)
    assert actual == [94, 86, 196, 0, 8, 0, 69, 127, 83, 52]


def test_sum():
    "Example to show sum of four tiles"
    assert 264 == sum([88, 96, 69, 11])


ipytest.run();
.........                                                                                    [100%]
9 passed in 0.01s

Genetic Programming Inspired Solution

The implementation creates a generation of tile layouts for each iteration until a layout is found that sums up to 264 across each row, column, and diagonal. There are few hyper parameters that were tweaked by hand while researching. The most important provided are what percentage of the top population to keep for mutations and how many new layouts to add. Finally, probabilities can be changed for the chances of flipping a tile in a mutation and if the tiles will be exchanged. Using 10% of the top tile layouts seemed to work the best, while generating a new population that is 90% of the original size by sampling with replacement mutating each. Finally, adding new scrambled layouts for the remaining 10%. Repeat until a solution is found.

It’s not perfect and it can get stuck in local minima. But, it is able find a solution in most cases in under 100,000 generations. Further research is needed for improvement.

g node0 Best 10% Rest 90% node1 Sample with Replacement and Exchange 90% New Population (Random Shuffle) 10% node0:f0->node1:f0 node2 Discard node0:f1->node2:f0
Figure 1: How New Generations Are Created.
import random
import math
from dataclasses import dataclass
from typing import Callable, Optional


def initial_pop(tiles: State = ALL_TILES, size: int = 1000) -> list[list[int]]:
    "Create a population of size of scrambled tiles"
    return [scramble(tiles) for _ in range(size)]


def sample_pop_best(pop: list[State], best: list[State]) -> list[State]:
    """
    Sample with replacement to create a new population of the best and add a population of new fully scrambled tiles.
    """
    return [
        exchange(each) for each in random.choices(best, k=len(pop) - len(best))
    ] + initial_pop(size=len(best))


@dataclass(frozen=True)
class Progress:
    """
    Simple holder for progress
    """

    generation: int
    mean_score: float
    best_score: float
    best_thus_far: State


def run(
    keep_percent: float = 0.1,
    progress_update: Callable[[Progress], None] = lambda progress: None,
    max_generations: int = 100_000,
) -> Optional[State]:
    """
    Simple solution using genetic programming concepts.
    """
    pop = initial_pop()
    generation = 1
    to_keep = int(len(pop) * keep_percent)
    while generation <= max_generations:
        scores = list(map(score, pop))
        top_best = sorted(zip(scores, pop), key=lambda each: each[0])[:to_keep]
        best_score = top_best[0][0]
        best = top_best[0][1]
        mean_score = sum(each[0] for each in top_best) / len(top_best)
        progress = Progress(generation, mean_score, mean_score, best)
        progress_update(progress)
        if math.isclose(best_score, 0.0):
            solution_index = scores.index(best_score)
            return pop[solution_index]
        pop = sample_pop_best(pop, [each[1] for each in top_best])
        generation += 1
    return None


@dataclass
class UpdatePrinter:
    """
    Simple object to be used for printing periodic progress
    """

    generation: int = 0
    how_often: int = 100

    def __call__(self, progress: Progress) -> None:
        self.generation = progress.generation
        if progress.generation % self.how_often == 0:
            print(
                "gen:",
                progress.generation,
                "mean:",
                round(progress.mean_score, 4),
                "best:",
                round(progress.best_score, 4),
                progress.best_thus_far,
            )

Sample run

%%time
seed = 1690419104
update = UpdatePrinter()
random.seed(seed)
solution = run(progress_update=update)
print("generations: ", update.generation)
print("solution: ", solution)
gen: 100 mean: 0.1423 best: 0.1423 [18, 89, 98, 61, 68, 91, 88, 16, 81, 19, 66, 99, 96, 69, 11, 86]
gen: 200 mean: 0.1404 best: 0.1404 [18, 91, 86, 69, 61, 88, 96, 18, 98, 11, 66, 89, 89, 66, 19, 91]
gen: 300 mean: 0.1322 best: 0.1322 [18, 89, 96, 61, 66, 88, 91, 19, 91, 18, 66, 86, 89, 69, 11, 98]
gen: 400 mean: 0.1441 best: 0.1441 [19, 89, 91, 66, 66, 88, 98, 11, 91, 18, 61, 89, 86, 69, 18, 96]
generations:  423
solution:  [18, 99, 86, 61, 66, 81, 98, 19, 91, 16, 69, 88, 89, 68, 11, 96]
CPU times: user 3.03 s, sys: 11.1 ms, total: 3.04 s
Wall time: 3.04 s

More solutions

ipytest.clean()

import pytest
import math
import random

SOLUTIONS = (
    (1669425992, [16, 98, 61, 89, 69, 81, 18, 96, 88, 66, 99, 11, 91, 19, 86, 68]),
    (1669426298, [68, 16, 81, 99, 89, 91, 66, 18, 96, 88, 19, 61, 11, 69, 98, 86]),
    (1690419104, [18, 99, 86, 61, 66, 81, 98, 19, 91, 16, 69, 88, 89, 68, 11, 96]),
)


@pytest.mark.parametrize("seed, expected", SOLUTIONS)
def test_solution(seed: int, expected: State):
    random.seed(seed)
    actual = run()
    assert actual == expected
    assert math.isclose(score(actual), 0.0, abs_tol=0.001)


ipytest.run();
...                                                                                          [100%]
3 passed in 109.71s (0:01:49)