Two dimensional Mastermind

Question

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.

Answer 1

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()[0]
          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()[0]
          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[3]),(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])

Answer 2

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[0] 

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[0] == 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[0] == 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[0]:
        for j in options[1]:
            for k in options[2]:
                for l in options[3]:
                    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[0] == 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):
                        new_column_possibilities[j].add(column_combo[j])
            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])