In a previous thread we made some progress fixing the one pre-existing backtrack algorithm, and did some basic testing to show that depth first searches can get results in practice. The purpose of this memo is to continue building towards a final implementation, which could end up as a deliverable to the Wolfram Function Repository.

Our previous backtrack algorithm assumed a level structure for the search space, where possible states on each level depend not on the current state of the search. This opens up a fine class of problems including some things about tilings and non-attacking queens; however, there is another obvious class of backtracking problems we would also like to solve some of. In the more general class, the set of states on each level does depend on the current state of the search (and likely also on history up to that point). In this case, we need the following alternative method for DepthFirstSearch:

DepthFirstIterate[iterator_, solutionQ_, count_
   ][state_] := Module[
  {nextStates = iterator[state],
   res = Nothing, allRes = {},
   ind = 1, len},
  len = Length[nextStates];
  If[len == 0,
   If[solutionQ[state],
    Return[{state}], Return[{Nothing}]]];
  While[ind <= len,
   Sow[DirectedEdge[state, nextStates[[ind]]]];
   res = DepthFirstIterate[iterator, solutionQ,
      count - Length[allRes]][
     nextStates[[ind]]];
   allRes = Join[allRes, res];
   If[Length[allRes] == count, Break[]];
   ind++]; allRes]


DepthFirstSearch[init_, iterator : (_Association | _Function),
  solutionQ_, count : (_Integer?Positive | All) : 1, False] := Block[
  {$RecursionLimit = Infinity},
  If[SameQ[Length[#], count, 1],
     First[#], #] &[
   DepthFirstIterate[
     iterator, solutionQ, count
     ][init]]]

A good test problem is to try and walk a chess knight around an $N \times N$ grid, visiting every square if possible, and sometimes even completing a Hamiltonian cycle. Per what is written on mathworld, this is not an easy task for a backtracking algorithm, so we will have to see what results we can get, if any.

First, let's define how a chess knight moves, and give an iterator for updating an $N \times N$ grid:

$knightMoves = Sort[Select[Tuples[{1, 2, -1, -2}, 2],
   UnsameQ @@ Abs[#] &]]

Out[]= {{-2, -1}, {-2, 1}, {-1, -2}, {-1, 2}, {1, -2}, {1, 
  2}, {2, -1}, {2, 1}}

IterateKnightWalk[board_] := With[{max = Max[board],
   pos = Position[board, Max[board]][[1]]},
  ReplacePart[board, # -> (max + 1)] & /@
   Select[pos + # & /@ RandomSample[$knightMoves],
    If[And @@ Thread[Dimensions[board] >= # > {0, 0}],
      board[[Sequence @@ #]] == 0, False] &]]

And if you look closely, you will see that the iterator contains some funny code. Usage of RandomSample allows that every time we run the backtracking search, potentially we will get a different answer (unless we set a SeedRandom). Furthermore, different answers take different amounts of time to compute, and (because the backtracker can easily get lost) similar answers are not guaranteed to compute in similar times:

data = Table[AbsoluteTiming[DepthFirstSearch[
      ReplacePart[ConstantArray[0, {8, 8}], {1, 1} -> 1],
      Function[IterateKnightWalk[#]], (True &),
      1] // Max], {5000}];

BoxWhiskerChart[Association[KeyValueMap[#1 -> #2[[All, 1]] &,
   KeySort@GroupBy[data, #[[2]] &]]],
 ChartLabels -> Keys[KeySort[GroupBy[data, #[[2]] &]]]]

It almost looks hopeful that the linear trend of minimal times could possibly continue to a complete board at 64, and in less than 0.01 seconds? Such a naive perspective ignores the multicomputational aspect of this plot, which we is plainly seen as times spiking up in the whiskers. Maximally slow times correspond to unlucky permutations via RandomChoice. To get a short path to a complete solution, we would have to get unbelievable lucky. In 5000 trials to 61 moves, we could only hit four times in less than 0.01 seconds:

Union[Table[TimeConstrained[AbsoluteTiming[DepthFirstSearch[
      ReplacePart[ConstantArray[0, {8, 8}], {1, 1} -> 1],
      Function[IterateKnightWalk[#]], (Max[#] == 62 &),
      1] // Max], .01], {5000}]] 


Out[] = {$Aborted, {0.006934`, 61}, {0.007279`, 61}, {0.007569`,   61}, {0.009356`, 61}}

The problem is pretty hard, but can we find anything on a smaller board? The $5 \times 5$ board is entirely computable in under 1 minute, and $6 \times 6$ is not too bad either. Hey look, I had some dumb luck on one of my searches:

AbsoluteTiming[DepthFirstSearch[
  ReplacePart[ConstantArray[0, {6, 6}], {3, 3} -> 1],
  Function[IterateKnightWalk[#]], (And[
     Max[#] == 36, MemberQ[$knightMoves,
      Position[#, 1][[1]] - Position[#, 36][[1]]]
     ] &), 1]]


{237.155305`, {{21, 6, 35, 2, 19, 4}, {36, 27, 20, 5, 34, 29}, {7, 22,
    1, 28, 3, 18}, {26, 15, 24, 11, 30, 33}, {23, 8, 13, 32, 17, 
   10}, {14, 25, 16, 9, 12, 31}}}

MatrixForm@{{21, 6, 35, 2, 19, 4}, {36, 27, 20, 5, 34, 29}, {7, 22, 1,
    28, 3, 18}, {26, 15, 24, 11, 30, 33}, {23, 8, 13, 32, 17, 
   10}, {14, 25, 16, 9, 12, 31}}

Not only is this a complete tour, it is also a Hamiltonian cycle on the underlying $6 \times 6$ knight graph because the 1 ends up a knight's move away from the 36. Just how lucky is it to find this result in only 4-5 minutes? I think that is a good type of question we should be trying to answer as part of the multicomputational paradigm. Anyone? For my part, I just tried to run the search again with TimeConstrained, and didn't find another such result. Moving on...

More manageable cases, where we can even make some interesting plots, occur when searching through the nodes of arbitrary Directed Acyclic Graphs. For example, let's define this as a search space:

ToDirectedAcyclicGraph = ResourceFunction["ToDirectedAcyclicGraph"];

SeedRandom[2414323];
With[{g0 =    WeaklyConnectedGraphComponents[RandomGraph[{300, 350}]][[1]]},
 g1 = Graph[ToDirectedAcyclicGraph[g0, VertexList[g0][[{1, 2}]]],
   VertexCoordinates -> Automatic, 
   GraphLayout -> "LayeredDigraphEmbedding",
   AspectRatio -> 1/2]] 

In this case we can iterate through levels just by following the edges:

edgeUpdate = Association[
   # -> VertexOutComponent[g1, #, {1}
       ] & /@ VertexList[g1]];

By default the search looks for a terminal node:

DepthFirstSearch[1, edgeUpdate, (True &)]
Out[]=217 

And we can see what this looks like by Reaping DirectedEdges:

HighlightGraph[g1,
 Graph[Flatten[Reap[
     DepthFirstSearch[1,
      edgeUpdate, (True &),
      1]][[2, 1]]]]]

And if we change the last argument to 2, 10 or All, we get:

The final graph should be the same as the entire out component of initial vertex 1. We can also notice that this graph has one deepest node at level 10, and plot the search tree for finding it:

HighlightGraph[g1,
 Graph[Flatten[Reap[
     DepthFirstSearch[1,
      edgeUpdate, (GraphDistance[g1, 1, #] == 10 &),
      1, TreeGraphQ -> True]][[2, 1]]]]]

That's neat, but can we get there from the other starting point? Just change the intial value and try to search again, the answer is yes:

HighlightGraph[g1,
 Graph[Flatten[Reap[
     DepthFirstSearch[23,
      edgeUpdate, (GraphDistance[g1, 1, #] == 10 &),
      1, TreeGraphQ -> True]][[2, 1]]]]]

That's it, basically. Because the algorithm can handle arbitrary DAGs and Iterators, it should be capable of solving a sufficiently wide range of backtracking problems. The last thing we have to worry about is whether a very deep graph could possibly hang a search. To avoid this we also want to introduce an optional cutoff for recursion depth, but this can probably wait until tomorrow, as well as writing the technical documentation...