# Performance tuning for game solving (peg solitaire | senku)

There's a game I saw at a friend's yesterday, that I often see at people's homes, but never for enough time to think on it too hard. It's called peg solitaire (thanks @R.M). So I came home and I wanted to find a solution in Mathematica, so I did the following

First, some visual functions. The game consists of a board with some slots that can either have a piece on it (black dot in this visual representation) or be empty (white dot)

empty=Circle[{0,0},0.3];
filled=Disk[{0, 0}, 0.3];

plotBoard[tab_]:=Graphics[GeometricTransformation[#1,TranslationTransform/@
Position[tab, #2]]&@@@{{empty, 0},{filled, 1}}, ImageSize->Small]


The starting board is the following.

tableroStart=({
{-1, -1, 1, 1, 1, -1, -1},
{-1, -1, 1, 1, 1, -1, -1},
{1, 1, 1, 1, 1, 1, 1},
{1, 1, 1, 0, 1, 1, 1},
{1, 1, 1, 1, 1, 1, 1},
{-1, -1, 1, 1, 1, -1, -1},
{-1, -1, 1, 1, 1, -1, -1}
});


-1 is used to represent places where there can't be any pieces. 0 for empty slots. 1 for slots with a piece on it.

So,

plotBoard[tableroStart] // Framed


Rules: Given a board such as the previous one, you can only move by "taking" a single piece, jumping over it. So, you take a piece, you choose one of the 4 straight directions, you jump over the adjacent piece and fall in an empty slot. The game is won by having only one last piece on the board. So, in the starting board, there are 4 possible moves, all symmetrical.

In this code, moves are represented by rules, so, {3, 4}->{3, 6} represents a move of the piece in coordinates {3, 4}, to coordinates {3, 6}, jumping over the piece at {3, 5} and taking it out of the board.

So, let's start programming.

This finds the possible moves towards some specified zero position

findMovesZero[tab_,pos_List]:=pos+#&/@(Join[#, Reverse/@#]&[Thread@{{0, 1, 3, 4}, 2}])//
Extract[ArrayPad[tab, 2],#]&//
Pick[{pos-{2, 0}, pos+{2, 0}, pos-{0, 2}, pos+{0, 2}},UnitStep[Total/@Partition[
#, 2]-2], 1]->pos&//Thread[#, List, 1]&


Lists all the possible moves given a board tab

i:findMoves[tab_]:=i=Flatten[#, 1]&[findMovesZero[tab, #]&/@Position[tab, 0]]


Given the board tab, makes the move

makeMove[tab_, posFrom_->posTo_]:=ReplacePart[tab , {posFrom->0, Mean[{posFrom, posTo}]->0,posTo->1}];


Now, the solving function

(* solve, given a board tab, returns a list of subsequent moves to win, or $Failed *) (* markTab is recursive. If a board is a success, marks it with$Success and makes all subsequent markTab calls return $NotNecessary *) (* If a board is not a success and doesn't have any more moves, returns$Failed. If it has moves, it just calls itself on every board,
saving the move made in the head of the new boards. I know, weird *)
Module[{$Success,$NotNecessary, parseSol, $guard, markTab}, markTab[tab_/;Count[tab, 1, {2}]===1]:=$Success/;!($guard=False)/;$guard;
i:markTab[tab_]:=With[{moves=findMoves[tab]},(i=If[moves==={}, $Failed,(#[markTab@makeMove[tab, #]]&/@moves)])]/;$guard;
markTab[tab_]/;!$guard:=$NotNecessary;

(* parseSol converts the tree returned by markTab into the list of moves until $Success, or in$Failed *)
parseSol[sol_]/;FreeQ[{sol}, $Success]:=$Failed;
parseSol[sol_]:=sol[[Apply[Sequence,#;;#&/@First@Position[sol, $Success]]]]//#/.r_Rule:>Null/;(Sow[r];False)&//Reap//#[[2, 1]]&; solve[tab_]:=Block[{$guard=True},parseSol@markTab@tab];
]


Solution visualization function

plotSolution[tablero_, moves_]:=
MapIndexed[Show[plotBoard[#1], Epilog->{Red,Dashed,Arrow[List@@First@moves[[#2]]]}]&, Rest@FoldList[makeMove[#, #2]&,tablero,moves]]//
Prepend[#, plotBoard[tablero]]&//Grid[Partition[#, 4, 4, 1, Null], Frame->All]&

(* Solves and plots *)
solveNplot = With[{sol=solve[#]},If[sol===$Failed,$Failed, plotSolution[#, sol]]]&;


In action:

solveNplot[( {
{-1, -1, 1, 1, 0, -1, -1},
{-1, -1, 1, 1, 1, -1, -1},
{1, 1, 0, 0, 0, 0, 0},
{1, 1, 0, 0, 0, 0, 0},
{1, 1, 0, 0, 0, 0, 0},
{-1, -1, 1, 1, 1, -1, -1},
{-1, -1, 1, 1, 1, -1, -1}
} )]


returns, after about a min's though,

So, the question is. How can we make it efficient enough so it can do the trick for an almost filled board like tableroStart?

The first move is actually always the same let alone symmetries so we could start a move ahead

• I love this game! It's called peg solitaire. I used to play it very often as a kid. The Wikipedia link also provides the shortest solution to the game. – rm -rf Jun 12 '12 at 5:39
• @R.M, SPOILER ALERT, I won't open the link! Shouldn't all solutions be equally short since you take no more and no less than one piece out every move? – Rojo Jun 12 '12 at 5:45
• Depends on your definition of 'move' and if you choose the wrong path in the decision tree. In the standard game, you can take several pieces in one move. For example, if you had a row 1 1 0 1 0 1 0, then you move the first element to the third position (at which point it can be captured by another piece from above or below). On the other hand, you could keep moving, resulting in 0 0 0 0 0 0 1 (and continue down/up and then back if the row position permits it). – rm -rf Jun 12 '12 at 6:14
• Shouldn't the final peg end up in the central hole? – Sjoerd C. de Vries Jun 12 '12 at 16:47
• @SjoerdC.deVries, I trust you more than my friend in this, so probably yes. This is just what my friend told me. It could be a nice follow-up – Rojo Jun 12 '12 at 17:32

## 2 Answers

### Preamble

Here is my first stab at it. This will not be the fastest possible solution (I hope to add some faster ones later), but even it will have no problems with your boards, including the full one you started with.

Before we dive into code, I will list the prerequisites for fast code in this case:

• Right choice of data structures
• Avoiding symbolic Mathematica overhead, which, sorry to say it, is just huge
• Avoiding copying in favor of direct modifications

### Code

Reproducing @Rojo's visualization functions to make this self-contained:

empty = Circle[{0, 0}, 0.3];
filled = Disk[{0, 0}, 0.3];

plotBoard[tab_] :=
Graphics[GeometricTransformation[#1,
TranslationTransform /@ Position[tab, #2]] & @@@
{{empty, 0}, {filled, 1}}, ImageSize -> Small]


I will start with your test board:

start =
{
{-1, -1, 1, 1, 0, -1, -1},
{-1, -1, 1, 1, 1, -1, -1},
{1, 1, 0, 0, 0, 0, 0},
{1, 1, 0, 0, 0, 0, 0},
{1, 1, 0, 0, 0, 0, 0},
{-1, -1, 1, 1, 1, -1, -1},
{-1, -1, 1, 1, 1, -1, -1}
}


First comes the optimized compiled function to find all possible steps for a given board:

getStepsC =
Compile[{{board, _Integer, 2}},
Module[{black = Table[{0, 0}, {Length[board]^2}], bctr = 0, i, j,
steps = Table[{{0, 0}, {0, 0}}, {Length[board]^2}], stepCtr = 0,
next, nnext
},
Do[
If[board[[i, j]] == 1, black[[++bctr]] = {i, j}],
{i, 1,Length[board]}, {j, 1, Length[board]}
];
black = Take[black, bctr];
Do[
Do[
next = pos + st;
nnext = pos + 2*st;
If[board[[next[[1]], next[[2]]]] == 1 &&
board[[nnext[[1]], nnext[[2]]]] == 0,
steps[[++stepCtr]] = {pos, nnext}
],
{st, {{1, 0}, {1, 1}, {0, 1}, {-1, 1},
{-1,0}, {-1, -1}, {0, -1}, {1, -1}}}
],
{pos, black}
];
Take[steps, stepCtr]],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
];


This function is expecting the board padded with -1-s, so that we don't have to check that the point belongs to the board. It will therefore also return cooridinates shifted by 1. It returns a list of sublists of starting and ending points for possible steps. Here is an example:

getStepsC[ArrayPad[start, 1, -1]]

  {{{2, 4}, {4, 4}}, {{2, 4}, {4, 6}}, {{2, 4}, {2, 6}}, {{2, 5}, {4, 5}},
{{2, 5}, {4, 7}}, {{4, 2}, {6, 4}}, {{4, 2}, {4, 4}}, {{5, 2}, {5, 4}},
{{6, 2}, {6, 4}}, {{6, 2}, {4, 4}}, {{8,4}, {6, 6}}, {{8, 4}, {6, 4}},
{{8, 5}, {6, 7}}, {{8, 5}, {6, 5}}, {{8, 6}, {6, 6}}, {{8, 6}, {6, 4}}}


Here is a function which helps to visualize all possible steps:

ClearAll[showPossibleSteps];
showPossibleSteps[brd_] :=
Show[plotBoard[brd],
Epilog ->
Map[{Red, Dashed, Arrow[# - {1, 1}]} &,
getStepsC[ArrayPad[brd, 1, -1]]]]


It pads the board with -1-s and subtracts 1 from both coordinates for the resulting steps. Using it, we get:

showPossibleSteps[start]


Next comes the main recursive function:

Clear[makeStep];
makeStep[steps : {step : {st_, end_}, prev_}, memoQ : (True | False) : False] :=
Module[{nblacks},
nblacks := Total@Clip[Flatten@board, {0, 1}];
If[nblacks == 1, Throw[steps, "Win"]];
If[memoQ && visited[board],
Return[]
];
board[[st[[1]], st[[2]]]] =
board[[(st[[1]] + end[[1]])/2, (st[[2]] + end[[2]])/2]] = 0;
board[[end[[1]], end[[2]]]] = 1;
If[nblacks == 1, Throw[steps, "Win"]];
Do[makeStep[{new, steps}, memoQ], {new, getStepsC[board]}];
If[memoQ, visited[board] = True];
board[[st[[1]], st[[2]]]] =
board[[(st[[1]] + end[[1]])/2, (st[[2]] + end[[2]])/2]] = 1;
board[[end[[1]], end[[2]]]] = 0;
];

makeStep[___] := Throw[\$Failed];


Few notes here: first, the board variable is not local to the body of makeStep (it is a global variable). Second, memoization can be switched on and off by the memoQ flag, and the related hash-table visited is also global. The above function is intended to be driven by the main one, not to be used in isolation. Last, note that the history of the previous steps is recorded in the linked list, which is an efficient way of doing this.

The way the function works is similarly to the @Rojo's code, but instead of collecting entire tree and then traversing it, it throws an exception at run-time as soon as the solution is found, and communicates the collected list of previous step via this exception. This allows the code to be memory-efficient.

Now, the main function:

Clear[getSolution];
getSolution[brd_, memoQ : (True | False) : False] :=
Block[{board = DeveloperToPackedArray@ArrayPad[brd, 1, -1], visited},
visited[_] = False;
Catch[
Do[makeStep[{new, {}}, memoQ], {new, getStepsC[board]}],
"Win"
]
];


Here are functions used for visualization:

ClearAll[showBoardStep];
showBoardStep[brd_, step_] :=
Show[plotBoard[brd], Epilog -> {Red, Dashed, Arrow[step]}];

ClearAll[toPlainListOfSteps];
toPlainListOfSteps[stepsLinkedList_] :=
Reverse@
Reap[
NestWhile[(Sow[First@# - {1, 1}]; Last[#]) &,
stepsLinkedList, # =!= {} &]
][[2, 1]];

ClearAll[showSolution];
showSolution[startBoard_, stepsLinkedList_] :=
Module[{b = startBoard},
Grid[Partition[#, 4, 4, 1, Null], Frame -> All] &@
MapAt[plotBoard, #, 1] &@
FoldList[
With[{st = #2[[1]], end = #2[[2]]},
b[[st[[1]], st[[2]]]] =
b[[(st[[1]] + end[[1]])/2, (st[[2]] + end[[2]])/2]] = 0;
b[[end[[1]], end[[2]]]] = 1;
showBoardStep[b, #2]] &,
b,
toPlainListOfSteps[stepsLinkedList]]];


What happens here is that I convert the linked list of steps to a plain list, and perform the relevant transformations on the board.

### Results and benchmarks

First, the test board, with and without the memoization:

getSolution[start]//Short//AbsoluteTiming

 {0.0585938,
{{{4,2},{4,4}},{{{4,5},{4,3}},{{{6,3},{4,5}},{{{7,4},{5,4}},
{{{8,6},{6,4}},<<1>>}}}}}
}

(stepList = getSolution[start,True])//Short//AbsoluteTiming

 {0.0419922,
{{{4,2},{4,4}},{{{4,5},{4,3}},{{{6,3},{4,5}},{{{7,4},{5,4}},
{{{8,6},{6,4}},<<1>>}}}}}
}


Note that the steps are reversed (last steps are shown first), and coordinates are shifted by 1. If you use

showSolution[start, stepList]


you get a sequence similar to what is displayed in the question.

Note that it only took a small fraction of a second to get the result (as opposed to a minute cited by @Rojo). Note also that memoization helped, but not dramatically so.

Now, the real deal:

(stepList0 = getSolution[tableroStart]);//AbsoluteTiming

 {18.7744141,Null}

(stepList = getSolution[tableroStart,True])//Short//AbsoluteTiming

 {2.0517578,{{{6,2},{6,4}},{{{6,5},{6,3}},{{{6,7},{6,5}},
{{{8,6},{6,6}},{{{8,4},{8,6}},<<1>>}}}}}}


Here memoization helps a great deal - we get an order of magnitude speedup. And here are the steps:

showSolution[tableroStart, stepList]


### Conclusions

This problem makes for a great case study, and is a very nice vehicle to study and analyze various performance issues as they reflect themselves in Mathematica. I have presented a straightforward (conceptually) implementation, whose main merit is not that the algorithm is particularly clever, but that it avoids some (but not all) serious performance pitfalls. Some other performance hits seem to be unavoidable, particularly those related to the top-level code being slow (makeStep function). This would have been different had Compile supported pass-by-reference and hash-tables (so that makeStep could be efficiently compiled).

As I said, this is not the fastest method, and I intend to add faster code later, but it illustrates main points. Note that the solution is essentially the same (conceptually) as what @Rojo did (except that I don't construct the full tree). What is really different is that frequent operations such as search for next steps are heavily optimized here (they take the most time), and also, I win big by mutating the board in place rather than copy it in the recursive invocations of makeStep. The result is a 3 orders of magnitude speed-up, and perhaps the solution has different computational complexity in general (although this is not yet clear to me).

Coming soon: Java port of this solution, prototyped entirely in Mathematica, which is another 20-30 times faster (according to my benchmarks).

• I don't understand this yet but +1 in confidence. – Mr.Wizard Jun 16 '12 at 15:45
• @Mr.Wizard Thanks :) – Leonid Shifrin Jun 16 '12 at 15:45
• +1. I'll dig into this soon, I was hoping you would get tempted by the performance tuning tag and come teach us, hehe. Anyway, it would probably be faster if you avoided the diagonal jumps, right? – Rojo Jun 16 '12 at 16:08
• @Rojo Somehow I wasn't reading your rules carefully enough, and assumed diagonal jumps are permitted. Which would be faster, hard to say, because diagonal jumps increase not only the number of all paths but also the number of valid solutions. But. generally, I think you are right. The modification to rule them out is trivial - one just has to remove half of the list of possible moves which in iterated over in the inner Do in the compiled function. – Leonid Shifrin Jun 16 '12 at 17:07
• @LeonidShifrin I am still surprised with how much faster your solution to the peg solitaire question is... I thought that given that mine also stops searching when it has found a solution, there wouldn't be a MMA solution 3 orders of magnitude faster! If that's all unnecessary symbolic overhead that I just threw out there without compunction, I'll have to start writing some mental notes – Rojo Jun 30 '12 at 12:52

Reader Beware: This does not solve the problem

I started trying to use graphs instead of lists with the idea of using isomorphisms to discard solution branches already explored.

I got a working solution, but too heavy to be used for the full problem. I was not able to find a way for efficiently identifying isomorphic states and so my solution is a recursive hog.

Anyway, as I am not going to spend more time with it for a while, I decided to post the code here, for the benefit of others trying the graph way.

Here it is:

(* manufacture vertex positions :) *)
places = MaxFilter[CrossMatrix[3], 1];
placesPos = Position[places, 1];
mPos = Max@placesPos;
pegNum = Length@placesPos;
g = Graph[Array[# &, pegNum], {}, VertexCoordinates -> placesPos, VertexLabels -> "Name", ImagePadding -> 10];

(*Allowable jumps *)

j1 = Select[Flatten[Table[Intersection[{{k, i}, {k, i + 1}, {k, i + 2}}, placesPos],
{k, mPos}, {i, mPos}], 1], Length@# == 3 &];
j2 = Select[Flatten[Table[Intersection[{{k, i}, {k + 1, i}, {k + 2, i}}, placesPos],
{k, mPos}, {i, mPos}], 1], Length@# == 3 &];
th = Thread[IntegerPart /@ PropertyValue[{g, #}, VertexCoordinates] & /@ VertexList[g] -> VertexList[g]];
j = Union[j1, j2] /. th;
jumps = UndirectedEdge @@@ (j[[All, {1, 3}]]);

(*test drive graph*)
g1 = Graph[Table[i, {i, pegNum}], jumps, VertexCoordinates -> placesPos,
VertexLabels -> "Name", ImagePadding -> 10];
(*jump over what vertex for each transition *)
(PropertyValue[{g1, #[[1]]}, "jOver"] = #[[2]]) & /@ Transpose[{jumps, j[[All, 2]]}];
(*Initially Occupied Vertex stock  Replace "7" with pegNum for a looong wait*)
(PropertyValue[{g1, #[[1]]}, "occ"] = #[[2]]) & /@ Table[{i, If[i <= 7, True, False]}, {i, pegNum}];

(*Utility funcs*)
SetAttributes[{freeV, jOver, pMoves, seeBoard, doMove, ret},  HoldFirst];
freeV[g_, x_] := ! PropertyValue[{g, x}, "occ"]; (*is the vertex free?*)
jOver[x_UndirectedEdge] := PropertyValue[{g1, x}, "jOver"]; (*Which vertex to jump over?*)
(*Select Possible moves at a certain graph state*)
pMoves[g_] := Select[jumps, ((freeV[g, #[[1]]]) != freeV[g, #[[2]]]) && (! freeV[g, jOver[#]]) &];
(*Utility for drawing  occupancy*)
seeBoard[g_] := Module[{},
vf[{xc_, yc_}, name_, {w_, h_}] :=If[freeV[g, name], {Blue, #}, {Red, #}] &@ Disk[{xc, yc}, Min@{w, h}];
Graph[Table[i, {i, pegNum}], jumps, VertexCoordinates -> placesPos,
VertexLabels -> "Name", ImagePadding -> 10, VertexShapeFunction -> vf, Frame -> True]];
(*perform a move> blank jOver vertex and traslate original*)
(*Note that we can't distinguish source & destination*)
(* Does not check if initial conditions are met*)
doMove[g_, x_UndirectedEdge] :=
(PropertyValue[{g, x[[1]]}, "occ"] = !PropertyValue[{g, x[[1]]}, "occ"];
PropertyValue[{g, x[[2]]}, "occ"] = !PropertyValue[{g, x[[2]]}, "occ"];
PropertyValue[{g, jOver[x]}, "occ"] = False;);
(*Test move*)
(*seeBoard[g1]
doMove[g1,1\[UndirectedEdge]9];*)
seeBoard[g1]
(*solving function*)
ret[g_, m_] := Module[{c := g},
If[(pMoves[g] != {}),
Module[{k = c}, (doMove[k, #]; ret[k, Append[m, #]])] & /@  pMoves[g]]; Sow[m]];
(*Check results. We dont distinguish between a->b and b->a yet!*)
l = (Reap@ret[g1, {}])[[2, 1]];
Length@l
TableForm@Select[l, (Length@# == (Max@(Length /@ l))) &]
`
• I thought about isomorphisms too. With the four-fold rotation symmetry of the starting position that would be a natural thing to do. I didn't go that way because it is my intuition that later on in the combinatorial explosion of this game symmetry will be less prominent and you don't gain sufficiently by pruning them. I may be wrong there. – Sjoerd C. de Vries Jun 13 '12 at 14:36
• @Sjoerd And there there are also reflection symmetries ... – Dr. belisarius Jun 13 '12 at 14:48
• @Sjoerd I think the Dihedral group D4/8 is a good candidate – Dr. belisarius Jun 13 '12 at 17:27
• There's also a deeper symmetry: the game is reversible and the roles of spaces and pieces can be reversed, too. – whuber Jun 14 '12 at 14:25