Bizarrely enough, I spent an undergrad summer in my Math degree studying this problem as well as coming up with algorithms to address it. First I will comment on the other answers.
Martin B: Correctly identified this as a Hamiltonian path problem. However, if the graph is a regular grid (as you two discussed in the comments) then a Hamiltonian path can trivially be found (snaking row-major order for example.)
agnorest: Correctly talked about the Hamiltonian path problem being in a class of hard problems. agnorest also alluded to possibly exploiting the graph structure, which funny enough, is trivial in the case of a regular grid.
I will now continue the discussion by elaborating on what I think you are trying to achieve. As you mentioned in the comments, it is trivial to find certain classes of "space-filling" non intersecting "walks" on a regular lattice/grid. However, it seems you are not satisfied only with these and would like to find an algorithm that finds "interesting" walks that cover your grid at random. But before I explore that possibility, I'd like to point out that the "non-intersecting" property of these walks is extremely important and what causes the difficulties of enumerating them.
In fact, the thing that I studied in my summer internship is called the Self Avoiding Walk. Surprisingly, the study of SAWs (self avoiding walks) is very important to a few sub-domains of modeling in physics (and was a key ingredient to the Manhattan project!) The algorithm you gave in your question is actually a variant on an algorithm known as the "growth" algorithm or Rosenbluth's algorithm (named after non other than Marshal Rosenbluth). More details on both the general version of this algorithm (used to estimate statistics on SAWs) as well as their relation to physics can be readily found in literature like this thesis.
SAWs in 2 dimensions are notoriously hard to study. Very few theoretical results are known about SAWs in 2 dimensions. Oddly enough, in higher than 3 dimensions you can say that most properties of SAWs are known theoretically, such as their growth constant. Suffice it to say, SAWs in 2 dimensions are very interesting things!
Reigning in this discussion to talk about your problem at hand, you are probably finding that your implementation of the growth algorithm gets "cut off" quite frequently and can not expand to fill the whole of your lattice. In this case, it is more appropriate to view your problem as a Hamiltonian Path problem. My approach to finding interesting Hamiltonian paths would be to formulate the problem as an Integer Linear Program, and add fix random edges to be used in the ILP. The fixing of random edges would give randomness to the generation process, and the ILP portion would efficiently compute whether certain configurations are feasible and if they are would return a solution.
The Code
The following code implements the Hamiltonian path or cycle problem for arbitrary initial conditions. It implements it on a regular 2D lattice with 4-connectivity. The formulation follows the standard sub-tour elimination ILP Lagrangian. The sub-tours are eliminated lazily, meaning there can be potentially many iterations required.
You could augment this to satisfy random or otherwise initial conditions that you deem "interesting" for your problem. If the initial conditions are infeasible it terminates early and prints this.
This code depends on NetworkX and PuLP.
"""
Hamiltonian path formulated as ILP, solved using PuLP, adapted from
https://projects.coin-or.org/PuLP/browser/trunk/examples/Sudoku1.py
Authors: ldog
"""
# Import PuLP modeler functions
from pulp import *
# Solve for Hamiltonian path or cycle
solve_type = 'cycle'
# Define grid size
N = 10
# If solving for path a start and end must be specified
if solve_type == 'path':
start_vertex = (0,0)
end_vertex = (5,5)
# Assuming 4-connectivity (up, down, left, right)
Edges = ["up", "down", "left", "right"]
Sequence = [ i for i in range(N) ]
# The Rows and Cols sequences follow this form, Vals will be which edge is used
Vals = Edges
Rows = Sequence
Cols = Sequence
# The prob variable is created to contain the problem data
prob = LpProblem("Hamiltonian Path Problem",LpMinimize)
# The problem variables are created
choices = LpVariable.dicts("Choice",(Vals,Rows,Cols),0,1,LpInteger)
# The arbitrary objective function is added
prob += 0, "Arbitrary Objective Function"
# A constraint ensuring that exactly two edges per node are used
# (a requirement for the solution to be a walk or path.)
for r in Rows:
for c in Cols:
if solve_type == 'cycle':
prob += lpSum([choices[v][r][c] for v in Vals ]) == 2, ""
elif solve_type == 'path':
if (r,c) == end_vertex or (r,c) == start_vertex:
prob += lpSum([choices[v][r][c] for v in Vals]) == 1, ""
else:
prob += lpSum([choices[v][r][c] for v in Vals]) == 2, ""
# A constraint ensuring that edges between adjacent nodes agree
for r in Rows:
for c in Cols:
#The up direction
if r > 0:
prob += choices["up"][r][c] == choices["down"][r-1][c],""
#The down direction
if r < N-1:
prob += choices["down"][r][c] == choices["up"][r+1][c],""
#The left direction
if c > 0:
prob += choices["left"][r][c] == choices["right"][r][c-1],""
#The right direction
if c < N-1:
prob += choices["right"][r][c] == choices["left"][r][c+1],""
# Ensure boundaries are not used
for c in Cols:
prob += choices["up"][0][c] == 0,""
prob += choices["down"][N-1][c] == 0,""
for r in Rows:
prob += choices["left"][r][0] == 0,""
prob += choices["right"][r][N-1] == 0,""
# Random conditions can be specified to give "interesting" paths or cycles
# that have this condition.
# In the simplest case, just specify one node with fixed edges used.
prob += choices["down"][2][2] == 1,""
prob += choices["up"][2][2] == 1,""
# Keep solving while the smallest cycle is not the whole thing
while True:
# The problem is solved using PuLP's choice of Solver
prob.solve()
# The status of the solution is printed to the screen
status = LpStatus[prob.status]
print "Status:", status
if status == 'Infeasible':
print 'The set of conditions imposed are impossible to solve for. Change the conditions.'
break
import networkx as nx
g = nx.Graph()
g.add_nodes_from([i for i in range(N*N)])
for r in Rows:
for c in Cols:
if value(choices['up'][r][c]) == 1:
nr = r - 1
nc = c
g.add_edge(r*N+c,nr*N+nc)
if value(choices['down'][r][c]) == 1:
nr = r + 1
nc = c
g.add_edge(r*N+c,nr*N+nc)
if value(choices['left'][r][c]) == 1:
nr = r
nc = c - 1
g.add_edge(r*N+c,nr*N+nc)
if value(choices['right'][r][c]) == 1:
nr = r
nc = c + 1
g.add_edge(r*N+c,nr*N+nc)
# Get connected components sorted by length
cc = sorted(nx.connected_components(g), key = len)
# For the shortest, add the remove cycle condition
def ngb(idx,v):
r = idx/N
c = idx%N
if v == 'up':
nr = r - 1
nc = c
if v == 'down':
nr = r + 1
nc = c
if v == 'left':
nr = r
nc = c - 1
if v == 'right':
nr = r
nc = c + 1
return nr*N+c
prob += lpSum([choices[v][idx/N][idx%N] for v in Vals for idx in cc[0] if ngb(idx,v) in cc[0] ]) \
<= 2*(len(cc[0]))-1, ""
# Pretty print the solution
if len(cc[0]) == N*N:
print ''
print '***************************'
print ' This is the final solution'
print '***************************'
for r in Rows:
s = ""
for c in Cols:
if value(choices['up'][r][c]) == 1:
s += " | "
else:
s += " "
print s
s = ""
for c in Cols:
if value(choices['left'][r][c]) == 1:
s += "-"
else:
s += " "
s += "X"
if value(choices['right'][r][c]) == 1:
s += "-"
else:
s += " "
print s
s = ""
for c in Cols:
if value(choices['down'][r][c]) == 1:
s += " | "
else:
s += " "
print s
if len(cc[0]) != N*N:
print 'Press key to continue to next iteration (which eliminates a suboptimal subtour) ...'
elif len(cc[0]) == N*N:
print 'Press key to terminate'
raw_input()
if len(cc[0]) == N*N:
break
