SET game odds simulation (MATLAB)

Question

I have recently found the great card came - SET. Briefly, there are 81 cards with the four features: symbol (oval, squiggle or diamond), color (red, purple or green), number (one, two or three) and shading (solid, striped or open). The task is to locate (from selected 12 cards) a SET of 3 cards, in which each of the four features is either all the same on each card or all different on each card (no 2+1 combination).

I've coded it in MATLAB to find a solution and to estimate odds of having a set in randomly selected cards.

Here is my code to estimate odds:

%% initialization
K = 12; % cards to draw
NF = 4; % number of features (usually 3 or 4)
setallcards = unique(nchoosek(repmat(1:3,1,NF),NF),'rows'); % all cards: rows - cards, columns - features
setallcomb = nchoosek(1:K,3); % index of all combinations of K cards by 3

%% test
tic
NIter=1e2; % number of test iterations
setexists = 0; % test results holder
% C = progress('init'); % if you have progress function from FileExchange
for d = 1:NIter
% C = progress(C,d/NIter);    

% cards for current test
setdrawncardidx = randi(size(setallcards,1),K,1);
setdrawncards = setallcards(setdrawncardidx,:);

% find all sets in current test iteration
for setcombidx = 1:size(setallcomb,1)
    setcomb = setdrawncards(setallcomb(setcombidx,:),:);
    if all(arrayfun(@(x) numel(unique(setcomb(:,x))), 1:NF)~=2) % test one combination
        setexists = setexists + 1;
        break % to find only the first set
    end
end
end
fprintf('Set:NoSet = %g:%g = %g:1\n', setexists, NIter-setexists, setexists/(NIter-setexists))
toc

100-1000 iterations are fast, but be careful with more. One million iterations takes about 15 hours on my home computer. Anyway, with 12 cards and 4 features I've got around 13:1 of having a set. This is actually a problem. The instruction book said this number should be 33:1. And it was recently confirmed by Peter Norvig. He provides the Python code, but I didn't test it yet.

So can you find an error? Any comments on performance improvement are welcome.

Answer 1

I tackled the problem writing my own implementation before looking at your code. My first attempt was very similar to what you already had :)

%# some parameters
NUM_ITER = 100000;  %# number of simulations to run
DRAW_SZ = 12;       %# number of cards we are dealing
SET_SZ = 3;         %# number of cards in a set
FEAT_NUM = 4;       %# number of features (symbol,color,number,shading)
FEAT_SZ = 3;        %# number of values per feature (eg: red/purple/green, ...)

%# cards features
features = {
    'oval' 'squiggle' 'diamond' ;    %# symbol
    'red' 'purple' 'green' ;         %# color
    'one' 'two' 'three' ;            %# number
    'solid' 'striped' 'open'         %# shading
};
fIdx = arrayfun(@(k) grp2idx(features(k,:)), 1:FEAT_NUM, 'UniformOutput',0);

%# list of all cards. Each card: [symbol,color,number,shading]
[W X Y Z] = ndgrid(fIdx{:});
cards = [W(:) X(:) Y(:) Z(:)];

%# all possible sets: choose 3 from 12
setsInd = nchoosek(1:DRAW_SZ,SET_SZ);

%# count number of valid sets in random draws of 12 cards
counterValidSet = 0;
for i=1:NUM_ITER
    %# pick 12 cards
    ord = randperm( size(cards,1) );
    cardsDrawn = cards(ord(1:DRAW_SZ),:);

    %# check for valid sets: features are all the same or all different
    for s=1:size(setsInd,1)
        %# set of 3 cards
        set = cardsDrawn(setsInd(s,:),:);

        %# check if set is valid
        count = arrayfun(@(k) numel(unique(set(:,k))), 1:FEAT_NUM);
        isValid = (count==1|count==3);

        %# increment counter
        if isValid
            counterValidSet = counterValidSet + 1;
            break           %# break early if found valid set among candidates
        end
    end
end

%# ratio of found-to-notfound
fprintf('Size=%d, Set=%d, NoSet=%d, Set:NoSet=%g\n', ...
    DRAW_SZ, counterValidSet, (NUM_ITER-counterValidSet), ...
    counterValidSet/(NUM_ITER-counterValidSet))

After using the Profiler to discover hot spots, some improvement can be made mainly by early-break'ing out of loops when possible. The main bottleneck is the call to the UNIQUE function. Those two lines above where we check for valid sets can be rewritten as:

%# check if set is valid
isValid = true;
for k=1:FEAT_NUM
    count = numel(unique(set(:,k)));
    if count~=1 && count~=3
        isValid = false;
        break   %# break early if one of the features doesnt meet conditions
    end
end

Unfortunately, the simulation is still slow for larger simulation. Thus my next solution is a vectorized version, where for each iteration, we build a single matrix of all possible sets of 3 cards from the hand of 12 drawn cards. For all these candidate sets, we use logical vectors to indicate what feature is present, thus avoiding the calls to UNIQUE/NUMEL (we want features all the same or all different on each card of the set).

I admit that the code is now less readable and harder to follow (thus I posted both versions for comparison). The reason being that I tried to optimize the code as much as possible, so that each iteration-loop is fully vectorized. Here is the final code:

%# some parameters
NUM_ITER = 100000;  %# number of simulations to run
DRAW_SZ = 12;       %# number of cards we are dealing
SET_SZ = 3;         %# number of cards in a set
FEAT_NUM = 4;       %# number of features (symbol,color,number,shading)
FEAT_SZ = 3;        %# number of values per feature (eg: red/purple/green, ...)

%# cards features
features = {
    'oval' 'squiggle' 'diamond' ;    %# symbol
    'red' 'purple' 'green' ;         %# color
    'one' 'two' 'three' ;            %# number
    'solid' 'striped' 'open'         %# shading
};
fIdx = arrayfun(@(k) grp2idx(features(k,:)), 1:FEAT_NUM, 'UniformOutput',0);

%# list of all cards. Each card: [symbol,color,number,shading]
[W X Y Z] = ndgrid(fIdx{:});
cards = [W(:) X(:) Y(:) Z(:)];

%# all possible sets: choose 3 from 12
setsInd = nchoosek(1:DRAW_SZ,SET_SZ);

%# optimizations: some calculations taken out of the loop
ss = setsInd(:);
set_sz2 = numel(ss)*FEAT_NUM/SET_SZ;
col = repmat(1:set_sz2,SET_SZ,1);
col = FEAT_SZ.*(col(:)-1);
M = false(FEAT_SZ,set_sz2);

%# progress indication
%#hWait = waitbar(0./NUM_ITER, 'Simulation...');

%# count number of valid sets in random draws of 12 cards
counterValidSet = 0;
for i=1:NUM_ITER
    %# update progress
    %#waitbar(i./NUM_ITER, hWait);

    %# pick 12 cards
    ord = randperm( size(cards,1) );
    cardsDrawn = cards(ord(1:DRAW_SZ),:);

    %# put all possible sets of 3 cards next to each other
    set = reshape(cardsDrawn(ss,:)',[],SET_SZ)';
    set = set(:);

    %# check for valid sets: features are all the same or all different
    M(:) = false;            %# if using PARFOR, it will complain about this
    M(set+col) = true;
    isValid = all(reshape(sum(M)~=2,FEAT_NUM,[]));

    %# increment counter if there is at least one valid set in all candidates
    if any(isValid)
        counterValidSet = counterValidSet + 1;
    end
end

%# ratio of found-to-notfound
fprintf('Size=%d, Set=%d, NoSet=%d, Set:NoSet=%g\n', ...
    DRAW_SZ, counterValidSet, (NUM_ITER-counterValidSet), ...
    counterValidSet/(NUM_ITER-counterValidSet))

%# close progress bar
%#close(hWait)

If you have the Parallel Processing Toolbox, you can easily replace the plain FOR-loop with a parallel PARFOR (you might want to move the initialization of the matrix M inside the loop again: replace M(:) = false; with M = false(FEAT_SZ,set_sz2);)

Here are some sample outputs of 50000 simulations (PARFOR used with a pool of 2 local instances):

» tic, SET_game2, toc
Size=12, Set=48376, NoSet=1624, Set:NoSet=29.7882
Elapsed time is 5.653933 seconds.

» tic, SET_game2, toc
Size=15, Set=49981, NoSet=19, Set:NoSet=2630.58
Elapsed time is 9.414917 seconds.

And with a million iterations (PARFOR for 12, no-PARFOR for 15):

» tic, SET_game2, toc
Size=12, Set=967516, NoSet=32484, Set:NoSet=29.7844
Elapsed time is 110.719903 seconds.

» tic, SET_game2, toc
Size=15, Set=999630, NoSet=370, Set:NoSet=2701.7
Elapsed time is 372.110412 seconds.

The odds ratio agree with the results reported by Peter Norvig.

Answer 2

Here's a vectorized version, where 1M hands can be calculated in about a minute. I got about 28:1 with it, so there might still be something a little off with finding 'all different' sets. My guess is that this is what your solution has trouble with, as well.

%# initialization
K = 12; %# cards to draw
NF = 4; %# number of features (this is hard-coded to 4)
nIter = 100000; %# number of iterations

%# each card has four features. This means that a card can be represented
%# by a coordinate in 4D space. A set is a full row, column, etc in 4D
%# space. We can even parallelize the iterations, at least as long as we
%# have RAM (each hand costs 81 bytes)
%# make card space - one dimension per feature, plus one for the iterations
cardSpace = false(3,3,3,3,nIter);

%# To draw cards, we put K trues into each cardSpace. I can't think of a
%# good, fast way to draw exactly K cards that doesn't involve calling
%# unique
for i=1:nIter
    shuffle = randperm(81) + (i-1) * 81;
    cardSpace(shuffle(1:K)) = true;
end

%# to test, all we have to do is check whether there is any row, column,
%# with all 1's
isEqual = squeeze(any(any(any(all(cardSpace,1),2),3),4) | ...
    any(any(any(all(cardSpace,2),1),3),4) | ...
    any(any(any(all(cardSpace,3),2),1),4) | ...
    any(any(any(all(cardSpace,4),2),3),1));
%# to get a set of 3 cards where all symbols are different, we require that
%# no 'sub-volume' is completely empty - there may be something wrong with this
%# but since my test looked ok, I'm not going to investigate on Friday night
isDifferent = squeeze(~any(all(all(all(~cardSpace,1),2),3),4) & ...
    ~any(all(all(all(~cardSpace,1),2),4),3) & ...
    ~any(all(all(all(~cardSpace,1),3),4),2) & ...
    ~any(all(all(all(~cardSpace,4),2),3),1));

isSet = isEqual | isDifferent;

%# find the odds
fprintf('odds are %5.2f:1\n',sum(isSet)/(nIter-sum(isSet)))

Answer 3

I found my error. Thanks Jonas for the hint with RANDPERM.

I used RANDI to randomly drawn K cards, but there is about 50% chance to get repeats even in 12 cards. When I substituted this line with randperm, I've got 33.8:1 with 10000 iterations, very close to the number in instruction book.

setdrawncardidx = randperm(81);
setdrawncardidx = setdrawncardidx(1:K);

Anyway, it would be interesting to see other approaches to the problem.