奇怪的是,我在数学学位的本科暑假里研究了这个问题,并提出了解决它的算法。首先,我将评论其他答案。
Martin B:正确地将其识别为哈密顿路径问题。但是,如果该图是一个规则网格(正如你们两个在评论中讨论的那样),那么可以很容易地找到一条哈密顿路径(例如,蛇行主要顺序。)
agnorest:正确地谈到哈密顿路径问题属于一类困难问题。agnorest 还暗示可能利用图结构,这很有趣,在规则网格的情况下是微不足道的。
我现在将通过详细说明我认为您正在努力实现的目标来继续讨论。正如您在评论中提到的,在常规格子/网格上找到某些类别的“空间填充”非相交“行走”是微不足道的。但是,您似乎不仅仅对这些感到满意,并且希望找到一种算法来找到随机覆盖您的网格的“有趣”步行。但在我探讨这种可能性之前,我想指出这些步行的“不相交”属性非常重要,是什么导致了枚举它们的困难。
事实上,我在暑期实习中学习的东西叫做自我避免步行。令人惊讶的是,SAW(自我避免行走)的研究对于物理学建模的一些子领域非常重要(并且是曼哈顿项目的关键成分!)您在问题中给出的算法实际上是算法的变体称为“增长”算法或 Rosenbluth 算法(以Marshal Rosenbluth命名)。有关该算法的一般版本(用于估计 SAW 的统计数据)以及它们与物理学的关系的更多详细信息,可以在本文这样的文献中轻松找到。
众所周知,二维的 SAW 很难研究。关于二维 SAW 的理论结果知之甚少。奇怪的是,在高于 3 维的情况下,您可以说 SAW 的大多数属性在理论上都是已知的,例如它们的生长常数。可以这么说,二维的 SAW 是非常有趣的东西!
在本次讨论中讨论您手头的问题时,您可能会发现增长算法的实现经常被“切断”,并且无法扩展以填充整个晶格。在这种情况下,将您的问题视为哈密顿路径问题更为合适。我寻找有趣的哈密顿路径的方法是将问题表述为整数线性程序,并添加要在 ILP 中使用的固定随机边。随机边缘的固定会给生成过程带来随机性,ILP 部分将有效地计算某些配置是否可行,如果可行,则返回解决方案。
编码
以下代码实现了任意初始条件的哈密顿路径或循环问题。它在具有 4 连通性的常规 2D 晶格上实现它。该公式遵循标准的子巡回消除 ILP 拉格朗日。子游览被延迟消除,这意味着可能需要多次迭代。
您可以增加它以满足您认为对您的问题“有趣”的随机或其他初始条件。如果初始条件不可行,它会提前终止并打印出来。
此代码取决于NetworkX和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