高效算法：基于评分的平面四色图着色

Question

我正在制作一个游戏，你需要使用四色定理为图像上色。相邻区域不能是同一种颜色。有四种颜色：红、绿、蓝、黄。每个红区域得10分，绿色得6分，蓝色得3分，黄色得1分。

我想要一个算法来计算任何给定图像的最高分数。我已经有了从图像中提取平面图的代码，它会给出每个区域的邻居列表。

到目前为止，我已经做了一个暴力实现，它检查所有可能的着色方案，但这样的复杂度是4的n次方，其中n是区域的数量。我可以采取的方法之一是尽可能优化这个搜索过程。

是否有更快的方法？我知道对于2种颜色，有一个线性时间的算法，但出于游戏设计原因，我通常不会生成可以用2种颜色涂色的图像。

谢谢 :)

编辑：如sascha的请求，这里有一些示例python字典，键是区域的id，列表是该区域的邻居列表

easy = {2: [4], 4: [2, 3, 14, 13], 3: [4], 14: [4], 13: [4]}

最高分数：46（我想）

（我的python暴力搜索0.1秒）

medium = {2: [4, 5, 6], 4: [2, 3], 3: [4, 18], 5: [2, 6], 6: [5, 2, 13, 18], 13: [6, 20, 21, 22], 18: [6, 3, 20, 22], 20: [18, 13], 22: [18, 13], 21: [13]}

最高分数：77

（我的python暴力搜索7.2秒）

hard = {2: [5, 6, 9], 5: [2, 4], 4: [5, 23], 6: [2, 7, 10], 3: [8, 16], 8: [3, 7, 12], 7: [6, 8, 10, 11], 9: [2, 10], 10: [6, 9, 7, 13, 14, 15, 17, 18], 11: [7, 12, 13], 12: [8, 11, 15, 16, 19], 13: [10, 11, 15], 14: [10, 15], 15: [10, 13, 12, 14, 17, 19], 16: [3, 12, 25, 27], 17: [10, 15, 18], 18: [10, 17, 19, 20], 19: [15, 18, 12, 27], 20: [18, 22, 24, 26, 27, 25], 22: [20], 23: [4, 24, 26], 24: [23, 20], 25: [16, 20], 26: [23, 20], 27: [19, 20, 16]}

(我的 Python 程序无法解决未知问题)

---

编辑2：

所以我完成了这个游戏，如果你感兴趣，可以在这里查看：here.

为了游戏的缘故，我意识到我只需要一个高分而不是绝对的最高分数（这是问题所要求的）。因此，我实现了贪心着色并运行了 10,000 次，每次洗牌图形并取得最好的得分结果。在所有小于 30 个区域的小板上，这将产生与暴力方法相同的结果，但它的时间复杂度更适用于较大的板。因此，它可能无法找到绝对最佳解决方案，但它总能找到一个非常好的解决方案。

非常感谢 @SaiBot 和 @sascha 的帮助 :)

- user4573386

---

获取所有可能的着色方案应该在O(2^n)内完成，因为在最坏情况下，你将有指数级别的可能性来给图像着色。不知道是否有帮助，但是在这篇论文中描述了一个在O(n^2)时间内对平面图进行四色着色的算法。也许可以对其进行调整，以根据您的评分模式优化着色方案。 - SaiBot

@sascha 谢谢你的建议，我已经添加了一些例子。如果你需要更多，只需提出 :) - user4573386

你有发布的图表的分数吗？比较一下会很好。 - SaiBot

@SaiBot，我在原帖中添加了一些内容，感谢您的帮助！ - user4573386

@SaiBot，我已经完成了这个游戏，如果你想试玩一下，我在原帖中（edit2）放了一个链接。非常感谢你的所有帮助。 - user4573386

Answer 1

这里介绍一种用Python进行简化整数规划的方法。

基本思路是利用现代高质量混合整数规划软件的惊人能力，而不必自己实现算法。 我们只需要定义模型（可能需要调整一些参数）！

请记住，（混合）整数规划通常是NP困难问题，并且我们假设这些启发式方法适用于我们的问题！

代码可能看起来有点丑陋，因为所使用的建模工具相当底层。 模型本身在结构上非常简单。

代码（python3; numpy、scipy、cylp + CoinOR Cbc求解器）

这是原型样式的代码，缺少最终解决方案的提取。由于这只是一个演示（您没有标记语言），这只是说明它可以是一个可行的方法。

from cylp.cy import CyClpSimplex
import itertools
import numpy as np
import scipy.sparse as sp
from timeit import default_timer as time

""" Instances """
# hard = {2: [4], 4: [2, 3, 14, 13], 3: [4], 14: [4], 13: [4]}
# hard = {2: [4, 5, 6], 4: [2, 3], 3: [4, 18], 5: [2, 6], 6: [5, 2, 13, 18], 13: [6, 20, 21, 22], 18: [6, 3, 20, 22], 20: [18, 13], 22: [18, 13], 21: [13]}
hard = {2: [5, 6, 9],
        5: [2, 4],
        4: [5, 23],
        6: [2, 7, 10],
        3: [8, 16],
        8: [3, 7, 12],
        7: [6, 8, 10, 11],
        9: [2, 10],
        10: [6, 9, 7, 13, 14, 15, 17, 18],
        11: [7, 12, 13],
        12: [8, 11, 15, 16, 19],
        13: [10, 11, 15],
        14: [10, 15],
        15: [10, 13, 12, 14, 17, 19],
        16: [3, 12, 25, 27],
        17: [10, 15, 18],
        18: [10, 17, 19, 20],
        19: [15, 18, 12, 27],
        20: [18, 22, 24, 26, 27, 25],
        22: [20],
        23: [4, 24, 26],
        24: [23, 20],
        25: [16, 20],
        26: [23, 20],
        27: [19, 20, 16]}

""" Preprocessing -> neighbor conflicts
    (remove dupes after sorting <-> symmetry

    Remark: for difficult use-cases one could try to think about special
    characteristics of the graph, like (not necessarily for this problem)
    chordal -> calc all max-cliques in P(olynomial-time) => pretty good convex-hull

    Here: just forbid conflicting-pairs (in each color-dimension).
"""

START_T = time()

conflicts = []
for key, vals in hard.items():
    for val in vals:
        conflicts.append((key, val))
conflicts_np = np.array(conflicts)
conflicts_np = np.sort(conflicts, axis=1)
conflicts_np = np.unique(conflicts_np, axis=0)

""" Preprocessing -> map IDs to gapless range [0-N)
"""
unique = np.unique(conflicts)
old2new = {}
new2old = {}
counter = itertools.count()
N = unique.shape[0]

for i in unique:
    new_id = next(counter)
    old2new[i] = new_id
    new2old[new_id] = i

conflicts_np = np.vectorize(old2new.get)(conflicts_np)

""" Sparse conflict matrix """
conflict_matrix = sp.coo_matrix((np.ones(conflicts_np.shape[0]*2),
                                (np.tile(np.arange(conflicts_np.shape[0]), 2),
                                 conflicts_np.ravel(order='F'))), shape=(conflicts_np.shape[0], N*4))
I, J, V = sp.find(conflict_matrix)

""" Integer Programming """

model = CyClpSimplex()
# 4 colors -> 4 binary vars per element in N
x = model.addVariable('x', N*4, isInt=True)
# scoring: linear-objective
model.objective = -np.hstack((np.full(N, 10), np.full(N, 6), np.full(N, 3), np.full(N, 1)))
# sub-opt way of forcing binary-constraints (from ints)
# (this awkward usage is due to problem with cylp in the past)
model += sp.eye(N*4) * x >= np.zeros(N*4)
model += sp.eye(N*4) * x <= np.ones(N*4)

# conflicts in each color-dimensions
# sub-opt numpy/scipy usage
for ind, i in enumerate(range(4)):
    if ind == 0:
        model += conflict_matrix * x <= 1
    else:
        shifted_conflicts = sp.coo_matrix((V,(I,J+(ind*N))), shape=(conflict_matrix.shape[0], N*4))
        model += shifted_conflicts * x <= 1

# force exactly one color per element
# sub-opt numpy/scipy usage
template = np.zeros(N*4)
template[0] = 1
template[N] = 1
template[2*N] = 1
template[3*N] = 1
all_color_dims = [sp.csc_matrix(np.roll(template, i).reshape(1,-1)) for i in range(N)]
model += sp.vstack(all_color_dims) *x == 1

cbcModel = model.getCbcModel() # Clp -> Cbc model / LP -> MIP
start_time = time()
status = cbcModel.solve()
end_time = time()

print("  CoinOR CBC used {:.{prec}f} secs".format(end_time - start_time, prec=3))
print("  Complete process used {:.{prec}f} secs".format(end_time - START_T, prec=3))

输出

Welcome to the CBC MILP Solver 
Version: 2.9.9 
Build Date: Jan 15 2018 

command line - ICbcModel -solve -quit (default strategy 1)
Continuous objective value is -200 - 0.00 seconds
Cgl0003I 0 fixed, 0 tightened bounds, 20 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 24 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 16 strengthened rows, 0 substitutions
Cgl0004I processed model has 153 rows, 100 columns (100 integer (100 of which binary)) and 380 elements
Cutoff increment increased from 1e-05 to 0.9999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of -194
Cbc0038I Before mini branch and bound, 100 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of -194 - took 0.00 seconds
Cbc0012I Integer solution of -194 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds)
Cbc0001I Search completed - best objective -194, took 0 iterations and 0 nodes (0.01 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from -194 to -194
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                -194.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.01

Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

  CoinOR CBC used 0.013 secs
  Complete process used 0.042 secs

结果

您的"large"实例在0.042秒内得到解决（包括子优化代码的完整时间），其中0.013秒用于核心求解器。当然，这只是一个实例，对其进行解释并不太科学！

结果与SaiBot的有趣定制解决方案（以及您的较小示例）相同！

（早期代码存在打分错误，因此我要求Saibot再次检查他的解决方案，现在我可以复现它！）

转移

MIP求解器应该在大多数架构和环境中都可用，甚至可能在移动设备上也可以使用（需要一些潜在的非平凡构建过程）。这些的建模/使用在某种程度上取决于建模系统和周围软件。

- sascha

---

谢谢@sascha，这真的很棒，速度非常快。 - user4573386

Answer 2

以下是我对该问题的尝试。我没有能够想出更好的时间复杂度，但优化了暴力算法。

我逐个处理节点，只允许着色，使得相邻节点不具有相同的颜色。

我为每个中间（非完整）着色添加了一个上限估计。为此，我假设每个未着色的节点将以最高得分的颜色进行着色（仅允许使用与已着色邻居不同的颜色）。因此，在上限计算中，两个尚未着色的相邻节点都可以被着色为“红色”。利用这个估计，我构建了一个分支界定算法，当上限估计仍低于当前最大值时，终止当前搜索路径。

小图的运行时间小于1毫秒，中等图为15毫秒，大图为3.2秒。三个图的结果分别为46、77和194。

import time
import copy

def upperBoundScore(graph, dfsGraphOder, dfsIndex, coloring, scoring, currentScore):
    maxAdditionalScore = 0;
    for i in range(dfsIndex, len(dfsGraphOder)):
        neighbourColors = {coloring[node] for node in graph[dfsGraphOder[i]]}
        possibleColors = {1, 2, 3, 4} - neighbourColors
        if len(possibleColors) < 1:  # if for one node no color is available stop
            return -1
        maxAdditionalScore += scoring[list(possibleColors)[0]]
    return currentScore+maxAdditionalScore

def colorRemainingGraph(graph, dfsGraphOder, dfsIndex, coloring, scoring, currentScore):
    global maxScore
    global bestColoring
    # whole graph colored
    if dfsIndex == len(dfsGraphOder):
        if currentScore > maxScore:
            maxScore = currentScore
            bestColoring = copy.deepcopy(coloring)
    # only proceed if current coloring can get better then best coloring
    elif upperBoundScore(graph, dfsGraphOder, dfsIndex, coloring, scoring, currentScore) > maxScore:
        neighbourColors ={coloring[node] for node in graph[dfsGraphOder[dfsIndex]]}
        possibleColors = list({1, 2, 3, 4} - neighbourColors)
        for c in possibleColors:
            coloring[dfsGraphOder[dfsIndex]] = c
            currentScore += scoring[c]
            colorRemainingGraph(graph, dfsGraphOder, dfsIndex+1, coloring, scoring, currentScore)
            currentScore -= scoring[c]
            coloring[dfsGraphOder[dfsIndex]] = 0

#graph = {2: [4], 4: [2, 3, 14, 13], 3: [4], 14: [4], 13: [4]}
#graph = {2: [4, 5, 6], 4: [2, 3], 3: [4, 18], 5: [2, 6], 6: [5, 2, 13, 18], 13: [6, 20, 21, 22], 18: [6, 3, 20, 22], 20: [18, 13], 22: [18, 13], 21: [13]}
graph = {2: [5, 6, 9], 5: [2, 4], 4: [5, 23], 6: [2, 7, 10], 3: [8, 16], 8: [3, 7, 12], 7: [6, 8, 10, 11], 9: [2, 10], 10: [6, 9, 7, 13, 14, 15, 17, 18], 11: [7, 12, 13], 12: [8, 11, 15, 16, 19], 13: [10, 11, 15], 14: [10, 15], 15: [10, 13, 12, 14, 17, 19], 16: [3, 12, 25, 27], 17: [10, 15, 18], 18: [10, 17, 19, 20], 19: [15, 18, 12, 27], 20: [18, 22, 24, 26, 27, 25], 22: [20], 23: [4, 24, 26], 24: [23, 20], 25: [16, 20], 26: [23, 20], 27: [19, 20, 16]}

# 0 = uncolored, 1 = red, 2 = green, 3 = blue, 4 = Yellow
scoring = {1:10, 2:6, 3:3, 4:1}
coloring = {node: 0 for node in graph.keys()}
nodeOrder = list(graph.keys())
maxScore = 0
bestColoring = {}

start = time.time()
colorRemainingGraph(graph, nodeOrder, 0, coloring, scoring, 0)
end = time.time()
print("Runtime: "+ str(end - start))
print("Max Score: "+str(maxScore))
print(bestColoring)

对于大图，得到的着色结果为(1=红色，2=绿色，3=蓝色，4=黄色):

{2: 1, 3: 1, 4: 1, 5: 2, 6: 2, 7: 1, 8: 2, 9: 2, 10: 3, 11: 2, 12: 1, 13: 1, 14: 1, 15: 4, 16: 2, 17: 2, 18: 1, 19: 2, 20: 2, 22: 1, 23: 2, 24: 1, 25: 1, 26: 1, 27: 1}

要验证算法输出的着色是否正确，可以使用以下代码，检查任何两个相邻节点是否具有相同的颜色。

def checkSolution(graph, coloring):
    validColoring=1
    for node in graph:
        for neighbour in graph[node]:
            if coloring[node] == coloring[neighbour]:
                print("wrong coloring found "+ str(node) + " and " + str(neighbour) + " have the same color")
                validColoring = 0
    if validColoring:
        print("Coloring is valid")

- SaiBot

---

这看起来是一个不错的答案。也许你想加入一些解决方案检查器。我的基于IP的求解器输出大数时是193，但我并不意味着这是正确的，也不意味着你的答案是错误的！ - sascha

@sascha 不错的想法！我添加了解决方案检查器以及大图的解决方案。如果有bug，也许你可以用它来找到它。输出的着色经过测试，所以194似乎是正确的答案。 - SaiBot

确实，你的分数看起来是正确的（现在可以复现）！ - sascha