# Two dimensional Mastermind

You have probably played the classic game of Mastermind with 4 pegs and 6 colours. It turns out that the codebreaker can always find the pattern in 5 moves or fewer.

Now consider the 2D version of the game. So we have a 4x4 grid of pegs and each peg can be one of 6 colours. Each guess involves making a 4x4 pattern of pegs. For each row and column, you will be told how many pegs have the correct colour and wrong position, and how many pegs have the correct colour and correct position. Obviously we can use the above classic algorithm to solve each row separately (in parallel), giving us a solution in 5 moves or fewer. However I am wondering whether we can reduce that to 4 moves or fewer? I believe it should be possible given that we now have more information coming from the columns.

I don't know the answer to this puzzle, so I am hoping that this amazing community can help me out.

• Note that the original version has 1296 possible codes, but this version has almost 3 trillion possible codes. The brute force approach in the Wikipedia article would probably still theoretically work, but nobody has the computing power necessary to make it feasible. Jul 21 at 19:30
• Since each guess gives out double the data points, three guesses doesn't seem to be out of the question by any means.
– Bass
Jul 21 at 21:28
• @Bass note that it's not as simple as "double the data" because the information is interlinked. For example a black hit in one column and one row pinpoints it exactly, while one in each of two rows and two columns narrows it down to two possibilities for which ones the black hits refer to. On the other hand a white hit in one column and one row means it is not the intersection point. Jul 21 at 22:03

What I'm going to describe is a rather gung-ho strategy which, frankly, shoudn't work, but to my surprise hasn't failed a single time in all of 100,000 trials I performed.

It finds the secret pattern after

3

or fewer trials. number of trials out of 100,000 where there was a given number of compatible patterns after two tests. number of annealing steps required for trials with high numbers of compatible patterns after two tests.

As pattern space is vast this doesn't rule out that there are worst case patterns defeating the strategy.

Strategy

The first two test patterns can---as found by trial and error---be chosen to be

     0023
1123
2355
2344

3411
3455
1102
5502


In particular, the second probing does not use any knowledge gained by the first.

The heuristics that went into the design of these patterns are that each pixel have 4 or 5 colours in the same row or column and access to one double colour in its row and in its column. Also patterns 1 and 2 should be as complementary as possible.

We can essentially by brute force find all patterns compatible with the scores (black and white pins) we received. Occasionally (~4%), there will be only one such pattern which then must be the solution, but typically it will be between a handful and few thousand (see chart) in which case we can tailor a third diagnostic pattern that separates all of them. I'm using a poor-man's simulated annealing (at least that's what I think I do) to find those third patterns and it works rather well.

Below is my Python code for those who would like to explore themselves:

import itertools as it
import numpy as np

# parameters

N,K = 6,4 # number of colours, number of spots
N_g = 100000 # total number of games to play
N_ann = 1000 # maximum number of annealing steps
N_rs = 1000 # maximum number of annealing restarts

# convention: we'll have to use K-dimensional and flat indexing side-by-side
# to minimise confusion we'll label _KD and _FL to indi

# precompute and tabulate stuff we will be need for gazillions of iterations
# (mostly ap_FL)

# for each pattern of K pins count how often each colour occurs
# (this is similar to numpy's bincount operation; hence the name bc_*)
# this is a K -> N dimensional map;
# therefore a full table will have dimensions N^{K+1}

bc_KD = np.zeros((K+1)*(N,),np.uint8)
for i in range(K):
# use einsum to access multidimensional diagonals
np.einsum("i...i->...i",bc_KD.swapaxes(0,i))[...] += 1
bc_FL = bc_KD.reshape(-1,N)

# for each pair of patterns of K pins compute "bp", the number of black pins
# this map is 2K -> 1 dimensional; full table size N^{2k}
bp_KD = np.zeros(2*K*(N,),np.uint8)
for i in range(K):
bp_KD0 = bp_KD.swapaxes(0,i)
bp_KD0 = bp_KD0.swapaxes(1,i+K)
np.einsum("ii...->i...",bp_KD0)[...] += 1
bp_FL = bp_KD.reshape(N**K,N**K)

# for each pair of patterns of K pins compute "tp", the total number of
# (black and white) pins
# this map is 2K -> 1 dimensional; full table size N^{2k}
tp_FL = np.minimum(bc_FL[:,None,:],bc_FL[None,:,:]).sum(-1)

# "to speed up later lookups encode white and black pin counts into a single
# integer "pp", "packed pins"
pp_FL = ((tp_FL*(tp_FL+1))>>1) + bp_FL
pp_KD = pp_FL.reshape(2*K*(N,))

# the total number of possible packed pin scores
# (note that for simplicity we do not use the fact that (K-1) x b + 1 x w
# is impossible)

# for a single row or column ...
n_pp = ((K+2)*(K+1))>>1
# ... for K rows or columns
N_pp = (n_pp)**K

# we'll also need some workspace
# this allocates 4 GB
ws = np.empty(1<<30,np.uint32)

del bc_KD,bc_FL,bp_KD,bp_FL

rng = np.random.default_rng()

# it is sometimes convenient to flatten N x N x ... x N (K factors) colour
# space coordinates into single integers. A K x K 2d pattern then still
# comprises K integers The following function computes the transpose.
def TKK(A):
return np.ravel_multi_index((*np.transpose(
np.unravel_index(A,K*(N,))),),K*(N,))

# given a list A of test patterns and the scores S (ST) for all rows (columns)
# compute all compatible patterns i.e. patterns that would yield the same
# scores when tested
# this is done by looking up all candidate rows and columns and then joining
# them discarding every combination that doesn't satisfy the constraint that
# row and column entries must agree at intersection points
# to minimise combinatorial explosion we alternate iterating over rows and
# columns, that way the constraint kicks in as early as possible
# this gives a dramatic speedup, well worth the overhead

# recursive function and ...
def _compat(CAR_KD,CAC_KD,selR_KD,selC_KD,n,m):
# iterate over nth row immediately filtering out all candidates that are
# not compatible with the columns already selected
for selR_KD[n] in CAR_KD[n][(CAR_KD[n][:,:m]==selC_KD[:m,n]).all(1)]:
if m == K:
yield selC_KD.copy()
else: # recurse, swapping the roles of rows and columns
yield from _compat(CAC_KD,CAR_KD,selC_KD,selR_KD,m,n+1)

# ... public entry point
def compat(A,SR,SC):
# allocate row and colmun slections
selR_KD,selC_KD = np.zeros((K,K),np.uint8),np.zeros((K,K),np.uint8)
AT = [TKK(a) for a in A]
AT = TKK(A)
# look up compatible rows and columns using the precomputed table
CAR_FL = [np.all([(pp_FL[a[i]] == s[i])
for a,s in zip(A,SR)],axis=0).nonzero()
for i in range(K)]
CAC_FL = [np.all([(pp_FL[a[i]] == s[i])
for a,s in zip(AT,SC)],axis=0).nonzero()
for i in range(K)]
CAR_KD = [np.transpose(np.unravel_index(c,K*(N,))) for c in CAR_FL]
CAC_KD = [np.transpose(np.unravel_index(c,K*(N,))) for c in CAC_FL]
yield from map(np.ravel_multi_index,_compat(
CAR_KD,CAC_KD,selR_KD,selC_KD,0,0),it.repeat(K*(N,)))

# given a list A of patterns and a test pattern B determine whether B separates
# A, i.e. whether each pattern in A scores differently against B

# this is the single most expensive operation
# to speed it up we use an arcane numpy feature
# when using advanced assignment with an index that contains duplicates while
# the data assigned do not then after reextraction using the same index
# the resulting vector will have the same number of elements as the original
# vectors but each group of elements that went through the same index value
# will all have been mapped to a single one of their original values the others
# being lost; which of the values is the one to be retained is undocumented

# for this trick to work we need an address space of size detemined by the
# largest possible index number
# in our case indices are full sets of scores; encoding each row and column
# independently this requires N_pp^2 x 2 x log2(N_pp) / 3 bytes of buffer

# we can, however make a small saving by using the fact that the total
# number of black pins must be the same over all rows and over all columns

# as an ugly but efficient optimisation allow transposes to be computed outside
# the function so it has to be done only once
def separates(B,BT,A,AT):
NA = len(A)
R = np.ravel_multi_index((*pp_FL[A,B].T,*pp_FL[AT,BT].T[:3],
tp_FL[AT,BT].T),(2*K-1)*(n_pp,) + (K+1,))
ws[R] = np.arange(NA)
coll, = (ws[R] != np.arange(NA)).nonzero()
return len(coll)

# this loop finds the third test pattern using biased random updates
def anneal(SR,SC,comp2):
comp2T = np.array(TKK(comp2))
# lowest number of collisions achieved so far
low = N_pp**2
# corresponding third test pattern
XR = np.zeros((K,K),np.uint8)
# annealing schedule; surely there is room for improvement
for t,T in enumerate(np.exp(-np.linspace(0,3,N_ann))):
# (candidate for) new third
NXR = XR.copy()
# randomize a random subset of pixels according to schedule
upd = rng.uniform(0,1,size=(K,K)) < T
NXR[upd] = (NXR[upd] + rng.integers(1,N,size=np.count_nonzero(upd),
dtype=np.uint8)) % N
# compress
NX = np.array(np.ravel_multi_index((*NXR.T,),K*(N,)))
NXT = np.array(np.ravel_multi_index((*NXR,),K*(N,)))
# check separation
code = separates(NX,NXT,comp2,comp2T)
if not code: # solution found
return True,t,XR
elif code < low: # improved but still insufficient third
low = code
XR = NXR
return False,t,XR

# hand crafted diagnostic patterns
AR = np.array([[[0,0,2,3],
[1,1,2,3],
[2,3,5,5],
[2,3,4,4]],
[[3,4,1,1],
[3,4,5,5],
[1,1,0,2],
[5,5,0,2]]],np.uint8)
# compress
A = np.ravel_multi_index((*AR.transpose(1,0,2),),K*(N,))
AT = np.ravel_multi_index((*AR.transpose(2,0,1),),K*(N,))
# place specific patterns you want to test in this list,
# they will be used before switching to random patterns
B_spec = np.empty((0,4,4),dtype=np.uint8)
# B_spec = [np.array([[4,1,5,0],
#                     [5,0,5,1],
#                     [5,4,1,4],
#                     [0,2,2,3]],dtype=np.uint8),
#           np.array([[5,0,4,0],
#                     [4,1,1,4],
#                     [2,4,3,1],
#                     [0,0,3,1]],dtype=np.uint8)]

# generate input patters
BR = rng.integers(0,N,size=(N_g,K,K),dtype=np.uint8)
BR[:len(B_spec)] = B_spec

# some stats
R_comp = [] # number of compatible patterns before third test
R_ann = []
R_rs = []
R_fail = 0 # number of games lost
fail = [] # collect more details about failures here
R_2 = 0
# main loop: every iteration is a game
for j,B,BT in zip(it.count(),
np.ravel_multi_index((*BR.transpose(1,0,2),),K*(N,)),
np.ravel_multi_index((*BR.transpose(2,0,1),),K*(N,))):
# compute scores against fixed test patterns
S = pp_FL[A,B]
SC = pp_FL[AT,BT]
# store list of compatible patterns
comp2 = np.array([*compat(A,S,SC)])
# retain some stats
R_comp.append(len(comp2))
if len(comp2) == 1:
print(f"({j}) fixed test patterns sufficient to identify {B}")
R_2 += 1
R_ann.append(0)
R_rs.append(0)
continue
for restart in range(N_rs):
code,t,X = anneal(S,SC,comp2)
if code:
R_ann.append(t)
R_rs.append(restart)
break
else: # finished schedule unsuccessfully
R_fail += 1
print(f"({j}) failed to separate {len(comp2)} patterns despite trying"
f" {N_rs} times using {N_ann} cooling steps")
# retain for post mortem
fail.append([B,comp2,X])

• I don't think you need to hide anything as a spoiler here. Also, what do you mean by saying some cases "appear already to be lost after the first two measurements." How many more tests do you need to do in order to solve them? The question is about how well it can do in the worst-case scenario, guaranteeing a win in x turns or fewer. It's great that this solves most cases in 3 steps, but if there are some cases where it will take you another 3 steps for a worst-case of 6, you're worse off than the naive solution of ignoring the columns. Jul 27 at 22:12
• I disagree. T he fact that it is so easy to reduce the expected number of steps below three says something about row-column synergy. By "lost" I meant because there are only finitely many black pin-white pin combinations (n_max, say) if before the third step we have more than n_max candidates left then we cannot finish with the third step if that were our ambition. So far I haven't encountered anything that would require more than four steps total, but that doesn't mean much. Jul 28 at 6:40
• I agree it's a good sign, but because the question is about what you can guarantee in a worst-case scenario it's weird that your analysis of the worst case is simply that it's "lost". I'd suggest updating your answer to include more info about what happens when you do need more steps. Jul 28 at 14:57
• With the new pair of fixed patterns and a more thorough annealing schedule (inlcuding restarts) I have yet to see a configuration that is not figured out in three. So far I have 1000 trials which admittedly is not all that many. I'll keep trying to speed up my code so I can start generating some proper statistics. I'll update the post once things settle down a bit. Jul 28 at 16:17
• @DmitryKamenetsky I've run 50,000 trials so far without a single failure. Maybe 5% don't require a third test. Typically number of patterns left after two tests are one or two digit numbers. The worst seen so far was between 2000 and 3000. It took more than 500 restarts to find a working third pattern there but eventually it was found. I'll add some charts to the post once I hit 100,000. Jul 29 at 9:36

I wrote a simple Python program that implements the classic solution, and the classic solution applied row-wise to the 2d version. Some things I've found:

• It takes a non-trivial amount of time to run. It took me a few minutes to run 100 random test cases. As I mentioned in a comment, there is no feasible way to fully extend the classic solution to 2d, as the computational requirements are ridiculously high.
• My implementation averages 4.42 moves to solve the single-row version.
• Solving the 2d version without taking columns into account usually takes 5 moves - it's about a 50% chance that an individual row takes 5 moves, so there's roughly a 1/16 chance that all 4 rows could be solved in fewer than 5 moves.

Next step: The easiest way to take columns into account is probably to just keep track of which colors are possible in which location. Then after each guess I update those possibilities. If any color is no longer possible in a particular location, any combinations that use that color in that location will be ruled out.

• 100 test runs resulted in 70 5's and 30 4's. The 30% chance of solving it in fewer than 5 tries is better than the naive 1/16, but worse than the classic solution with a single row.

At this point I can't say if it's possible to guarantee a solution in 4 or fewer. I've tried a couple initial guesses, but none seem to do better than the one I've got in the code below. Note that this says more about my strategy than it does about the problem itself - when trying to apply the classic solution to the 2d version it is nontrivial to also apply information gained from the columns.

Here's my code:

from collections import Counter
from random import randint

def score_row(guess, actual):
black_hits = 0
white_hits = 0

guess_counter = Counter(guess)
actual_counter = Counter(actual)
for color, count in guess_counter.items():
white_hits += min(count, actual_counter[color])

for i, j in zip(guess, actual):
if i == j:
black_hits += 1
white_hits -= 1

return (black_hits, white_hits)

def combo_iterator():
# generate all possible combos (0,0,0,0) through (5,5,5,5)
for i in range(6):
for j in range(6):
for k in range(6):
for l in range(6):
yield (i, j, k, l)

def determine_next_guess(guess_choices, possible_codes):
for guess_choice in guess_choices:
if guess_choice in possible_codes:
return guess_choice
return guess_choices

def run_normal_game():
actual_code = tuple(randint(0,5) for i in range(4))
turns = 0
next_guess = (0, 0, 1, 1)
possible_codes = set(combo_iterator())

while True:
turns += 1
# Play guess to get response of pegs
guess_score = score_row(next_guess, actual_code)
if guess_score == 4:
return turns

# Remove from S any codes that don't match response
codes_to_remove = []
for code in possible_codes:
if score_row(next_guess, code) != guess_score:
codes_to_remove.append(code)
for code in codes_to_remove:
possible_codes.remove(code)

# Determine next guess
max_minimum = 0
guess_choices = []
for guess_choice in combo_iterator():
hit_counts = Counter()
for possible_code in possible_codes:
hit_counts.update((score_row(guess_choice, possible_code),))
current_score = len(possible_codes) - max(hit_counts.values())
if current_score >= max_minimum:
if current_score > max_minimum:
max_minimum = current_score
guess_choices = []
guess_choices.append(guess_choice)
next_guess = determine_next_guess(guess_choices, possible_codes)

def run_2d_game():
actual_code = tuple(tuple(randint(0,5) for i in range(4)) for j in range(4))
turns = 0
next_guess = [(0, 0, 1, 1), (0, 0, 1, 1), (1, 1, 2, 2), (1, 1, 2, 2)]
possible_codes = tuple(set(combo_iterator()) for i in range(4))

while True:
turns += 1
# Play guess to get response of pegs
guess_scores = [score_row(next_guess[i], actual_code[i]) for i in range(4)]
if all(guess_score == 4 for guess_score in guess_scores):
return turns

for i in range(4):
# Remove from S any codes that don't match response
codes_to_remove = []
for code in possible_codes[i]:
if score_row(next_guess[i], code) != guess_scores[i]:
codes_to_remove.append(code)
for code in codes_to_remove:
possible_codes[i].remove(code)

# Determine next guess
max_minimum = 0
guess_choices = []
for guess_choice in combo_iterator():
hit_counts = Counter()
for possible_code in possible_codes[i]:
hit_counts.update((score_row(guess_choice, possible_code),))
current_score = len(possible_codes[i]) - max(hit_counts.values())
if current_score >= max_minimum:
if current_score > max_minimum:
max_minimum = current_score
guess_choices = []
guess_choices.append(guess_choice)
next_guess[i] = determine_next_guess(guess_choices, possible_codes[i])


Here's the 2nd part (taking columns into account):

def column_combo_iterator(options):
for i in options:
for j in options:
for k in options:
for l in options:
yield (i, j, k, l)

def run_2d_game_with_columns():
actual_code = tuple(tuple(randint(0,5) for i in range(4)) for j in range(4))
turns = 0
next_guess = [(0, 0, 1, 1), (0, 0, 1, 1), (1, 1, 2, 2), (1, 1, 2, 2)]
possible_codes = tuple(set(combo_iterator()) for i in range(4))
column_info = tuple(tuple(set(range(6)) for i in range(4)) for j in range(4))

while True:
turns += 1
# Play guess to get response of pegs
guess_scores = [score_row(next_guess[i], actual_code[i]) for i in range(4)]
if all(guess_score == 4 for guess_score in guess_scores):
return turns

for i in range(4):
new_column_possibilities = [set() for j in range(4)]
guessed_column = tuple(next_guess[j][i] for j in range(4))
column_score = score_row(guessed_column, tuple(actual_code[j][i] for j in range(4)))
for column_combo in column_combo_iterator(column_info[i]):
if score_row(guessed_column, column_combo) == column_score:
for j in range(4):
for j in range(4):
column_info[i][j].intersection_update(new_column_possibilities[j])

for i in range(4):
# Remove from S any codes that don't match response
codes_to_remove = []
for code in possible_codes[i]:
if score_row(next_guess[i], code) != guess_scores[i] or any(code[j] not in column_info[j][i] for j in range(4)):
codes_to_remove.append(code)
for code in codes_to_remove:
possible_codes[i].remove(code)

# Determine next guess
max_minimum = 0
guess_choices = []
for guess_choice in combo_iterator():
hit_counts = Counter()
for possible_code in possible_codes[i]:
hit_counts.update((score_row(guess_choice, possible_code),))
current_score = len(possible_codes[i]) - max(hit_counts.values())
if current_score >= max_minimum:
if current_score > max_minimum:
max_minimum = current_score
guess_choices = []
guess_choices.append(guess_choice)
next_guess[i] = determine_next_guess(guess_choices, possible_codes[i])

• Thank you for your work. This is very interesting. It is nice that you managed to add the columns information into the classic algorithm. It seems that it only makes a marginal improvement though. Jul 29 at 0:00