1. 项目概述与问题拆解最近在整理一些经典的逻辑与编程练习题时又遇到了这个“20棵树每行4棵最多能排多少行”的问题。这听起来像是个简单的排列组合或者小学数学题但实际上它背后牵扯到的是计算几何中一个非常有趣且历史悠久的话题——果园问题Orchard-planting problem或者更具体地说是“西尔维斯特直线问题”的一个变种。很多朋友第一次接触时可能会尝试手动画图或者用排列组合公式去计算但很快就会陷入困境因为点的位置是连续的可能的排列方式是无限的。这正是计算机模拟可以大显身手的地方我们不需要去证明一个完美的数学解而是可以通过程序去探索、去逼近在有限的尝试中找到当前最优的布局。今天我就手把手地带大家用Python从零开始构建一个模拟程序来探索这个问题的答案并深入理解其中的算法思想和优化技巧。简单来说我们的目标是在平面上放置20个点代表树寻找一种布局使得至少包含3个点的直线数量尽可能多并且每条这样的直线上恰好有4个点。这里“至少3个点”是构成直线的条件而我们的目标是每条线4个点所以实际上我们寻找的是所有恰好包含4个点的直线。我们将编写一个程序能够自动生成点的布局、检测所有可能的4点共线情况并通过迭代优化来尝试增加共线直线的数量。这个过程不仅锻炼编程能力更能让我们直观地感受组合优化问题的魅力。无论你是Python新手想找一个有挑战性的练手项目还是对算法优化感兴趣的老手相信都能从中获得启发。2. 核心思路与算法设计面对这样一个问题直接进行数学求解是极其困难的。因此我们转向计算模拟和启发式搜索。核心思路可以概括为生成候选点集 - 检测所有4点共线组合 - 评估目标函数共线直线数 - 迭代优化点集以提升目标函数。2.1 问题建模与关键定义首先我们需要将问题转化为计算机可以处理的形式。点用一个二维坐标(x, y)表示一棵树的位置。为了简化我们可以将点限制在一个有限区域内比如一个100x100的整数坐标网格上或者使用浮点数在单位正方形[0,1] x [0,1]内随机生成。直线与共线判定如何判断四个点是否在同一直线上这是算法的核心。在离散网格上我们可以计算斜率但需要处理垂直线斜率无穷大和浮点数精度问题。更稳健的方法是使用向量共线性原理。对于三个点A, B, C可以通过计算向量AB和AC的叉积在二维中体现为行列式是否为零来判断它们是否共线。对于四个点我们可以先判断前三个点共线再判断第四个点是否也在这条线上即满足同样的直线方程。注意由于计算机使用浮点数直接判断行列式 0是不可靠的。我们必须引入一个容差tolerance例如abs(行列式值) 1e-10来判断是否“近似共线”。这个容差的选择需要谨慎太小会漏判太大会误判。目标函数对于一个包含20个点的集合S定义其目标函数值f(S)为从S中找出所有满足“恰好包含4个点且这4点共线”的子集这些子集所对应的不同直线的数量。注意同一条直线上可能有超过4个点比如5个但我们只关心恰好4个点构成的行。如果一条线上有5个点那么从中任取4个点都共线但这只算作一条有效的行而不是C(5,4)5条。因此我们需要对检测到的4点组合进行去重识别出它们所属的唯一直线。2.2 算法框架选择我们不可能遍历所有可能的20个点的布局那是无限种。因此必须采用启发式搜索算法。这里我选择一种简单但有效的迭代优化框架结合局部搜索和随机重启。初始化随机生成20个点的初始布局例如在单位正方形内随机均匀分布。评估编写函数evaluate(points)输入20个点的列表返回检测到的有效行4点共线且唯一的数量以及所有检测到的4点组合用于后续分析。邻居生成局部搜索定义当前布局的一个“邻居”。一个简单的邻居生成方式是随机选择一个点在其当前位置附近的一个小邻域内如半径为delta的圆内随机扰动其位置得到一个新的布局。搜索策略采用模拟退火Simulated Annealing或爬山法Hill Climbing的变种。爬山法生成一个邻居布局计算其目标函数值。如果新值优于当前值则接受这个邻居作为新的当前布局否则拒绝。重复此过程直到连续多次迭代没有改进。模拟退火在爬山法的基础上允许以一定的概率接受比当前解差的邻居这个概率随着“温度”的降低而减小。这有助于跳出局部最优解。对于本问题模拟退火通常能取得更好的效果。随机重启由于搜索空间复杂算法很容易陷入局部最优。因此在单次搜索收敛后我们保留找到的最好解然后完全随机生成一个新的初始布局重新开始搜索。重复这个过程多次最后输出所有随机重启中找到的全局最优解。2.3 共线检测算法优化暴力检测所有4点组合C(20,4)4845种并判断共线在每次评估时都进行计算量是可以接受的。但我们可以进一步优化提前终止对于三个点如果它们不共线那么包含这三个点的任何四个点集合也一定不共线对于第四个点的任意选择。但这在暴力枚举中优化效果有限。基于斜率和截距的哈希判断四个点(p1, p2, p3, p4)共线等价于它们满足同一条直线方程y kx b或x c。我们可以先计算每两个点确定的直线用标准化后的(k, b)或(c)表示然后检查是否有四个点共享同一条直线。具体步骤遍历所有点对(i, j)计算它们确定的直线L。用一个字典哈希表记录每条直线L上包含哪些点。遍历结束后检查字典中哪些直线对应的点集大小 4。对于点集大小 4的直线其包含的4点组合数非常多但我们只关心这条直线本身。所以每一条这样的直线只要其包含的点数4就为最终行数贡献1。 这种方法比暴力枚举4点组合更高效尤其是当点数增多时。复杂度从 O(n^4) 降到了 O(n^2)。但需要注意浮点数精度问题直线参数(k, b)需要量化如四舍五入到小数点后若干位才能作为字典的键。在本教程中为了代码清晰和易于理解我们先实现暴力枚举法。在后续优化部分我们再探讨哈希法的实现。3. 基础版本实现暴力枚举与评估让我们开始动手编码。首先我们实现最核心的共线检测和目标函数评估。3.1 环境准备与点表示我们将使用Python的标准库math、random和itertools。不需要安装额外包非常适合新手。import random import math import itertools from typing import List, Tuple, Set Point Tuple[float, float] # 类型别名表示一个点 (x, y)我们定义一个函数来生成随机初始点集。这里选择在[0, 100] x [0, 100]的正方形区域内生成整数坐标方便观察和调试。你也可以使用random.uniform(0, 100)生成浮点数。def generate_random_points(num_points: int 20, coord_range: Tuple[int, int] (0, 100)) - List[Point]: 生成指定数量的随机点。 Args: num_points: 点的数量默认为20。 coord_range: 坐标范围默认为(0, 100)。 Returns: 一个包含num_points个随机点的列表。 low, high coord_range return [(random.randint(low, high), random.randint(low, high)) for _ in range(num_points)]3.2 共线判断函数实现一个稳健的共线判断函数用于判断三个点是否共线以及第四个点是否在前三个点确定的直线上。def are_collinear(p1: Point, p2: Point, p3: Point, tolerance: float 1e-10) - bool: 判断三个点p1, p2, p3是否共线。 使用三角形面积公式行列式|x1(y2-y3) x2(y3-y1) x3(y1-y2)| / 2。 面积为零则共线。 Args: p1, p2, p3: 三个点。 tolerance: 容差用于处理浮点数精度误差。 Returns: True如果三点共线在容差范围内。 x1, y1 p1 x2, y2 p2 x3, y3 p3 # 计算行列式值面积的2倍 area2 abs(x1 * (y2 - y3) x2 * (y3 - y1) x3 * (y1 - y2)) return area2 tolerance def is_point_on_line(p1: Point, p2: Point, px: Point, tolerance: float 1e-10) - bool: 判断点px是否在由点p1和p2确定的直线上。 原理向量(p1-px)与向量(p1-p2)共线。 等价于判断三个点p1, p2, px共线。 return are_collinear(p1, p2, px, tolerance)3.3 目标函数评估暴力枚举法现在实现评估函数。暴力枚举所有4点组合检查它们是否共线。然后我们需要对共线的4点组合进行去重得到唯一的直线。去重是关键。两条4点共线组合可能代表同一条直线吗是的例如直线L上有5个点A,B,C,D,E。那么组合 {A,B,C,D}、{A,B,C,E}、{A,B,D,E} 等都代表同一条直线L。我们需要将它们合并。一个简单但有效的去重方法是对于每一个检测到的4点共线组合我们计算这条直线的“签名”signature。我们可以用直线上任意两点计算斜率和截距但垂直线需要特殊处理。更稳健的“签名”是将这条直线上的所有点当前是4个按坐标排序取最远的两个点或者首尾两点来确定直线然后对这个确定的直线进行标准化表示。这里我们采用一个更直接的方法点集排序法。当我们找到一个4点共线组合时我们得到的是4个点。我们把这4个点按照x坐标如果x相同则按y坐标排序。用排序后的点列表作为该直线的“代表”。但注意如果一条直线包含5个点我们可能会得到多个不同的4点子集它们的排序列表不同但属于同一直线。因此更好的方法是当我们找到一条包含4个点的直线时我们找出这条直线上当前所有的点。对于暴力枚举我们暂时只知道这4个点。我们可以先记录这4个点。在评估函数的最后我们需要合并那些代表同一直线的点集。判断两个点集是否在同一直线上可以检查它们是否共享至少两个公共点并且所有点都共线。如果是则合并两个点集。合并后我们得到若干个点集每个点集代表一条直线并且点集的大小4。那么最终的行数就是这些点集的数量。根据这个思路我们实现第一版的评估函数def evaluate_bruteforce(points: List[Point], collinear_tol: float 1e-10) - Tuple[int, List[Set[int]]]: 暴力评估函数计算给定点集中有多少条不同的直线恰好包含至少4个点我们目标是4个。 注意这里返回的是包含4个点的直线数量。我们最终目标是最大化这个数。 Args: points: 点列表。 collinear_tol: 共线判断容差。 Returns: (行数, 直线列表) 直线列表中的每个元素是一个点索引的集合代表一条直线上的点。 n len(points) all_lines [] # 存储找到的直线每条直线用一个点索引集合表示 # 枚举所有4点组合 for quad in itertools.combinations(range(n), 4): i, j, k, l quad p1, p2, p3, p4 points[i], points[j], points[k], points[l] # 判断四点是否共线先判断前三点再判断第四点 if are_collinear(p1, p2, p3, collinear_tol) and is_point_on_line(p1, p2, p4, collinear_tol): # 找到一个4点共线组合 current_line_indices set(quad) # {i, j, k, l} # 检查这个组合是否与已找到的某条直线是同一根即直线上有更多点 merged False for line in all_lines: # 如果当前组合与已有直线共享至少2个点则认为是同一直线合并 if len(current_line_indices line) 2: # 还需要验证合并后的所有点是否依然共线理论上应该但浮点容差下需确认 # 简单起见假设共享2点即共线直接合并 line.update(current_line_indices) merged True break if not merged: all_lines.append(current_line_indices) # 合并后我们可能得到一些点集大小4的直线。我们只关心直线本身。 # 过滤掉点数小于4的直线理论上不会出现但安全起见 valid_lines [line for line in all_lines if len(line) 4] num_rows len(valid_lines) return num_rows, valid_lines这个函数能工作但效率不高且合并逻辑在浮点容差和复杂情况下可能不可靠。但它为我们提供了一个起点和验证基准。3.4 简单测试与可视化让我们写一个简单的测试并尝试用字符画来可视化结果对于简单布局可行。def print_points_and_lines(points: List[Point], lines: List[Set[int]], coord_range(0,100)): 简陋的文本可视化将坐标区域网格化打印点用数字索引和直线用字母表示。 仅适用于小范围整数坐标和少量点用于调试。 low, high coord_range size high - low 1 # 创建网格 grid [[. for _ in range(size)] for __ in range(size)] # 绘制点 for idx, (x, y) in enumerate(points): ix, iy int(round(x)) - low, int(round(y)) - low if 0 ix size and 0 iy size: # 用索引的个位数表示点超过9用字母 ch str(idx) if idx 10 else chr(ord(A) idx - 10) grid[iy][ix] ch # 打印网格 (y轴从上到下递减符合常规坐标系) print(Points layout (indices):) for row in reversed(grid): # 反转Y轴以便从下往上打印 print( .join(row)) print(f\nFound {len(lines)} lines with 4 points:) for i, line_indices in enumerate(lines): points_on_line [points[idx] for idx in line_indices] print(f Line {i}: indices {sorted(line_indices)}, points {points_on_line}) # 简单测试 if __name__ __main__: random.seed(42) # 固定随机种子使结果可复现 test_points generate_random_points(10, (0, 10)) # 先用10个点在小范围测试 print(Generated points:, test_points) num_rows, lines evaluate_bruteforce(test_points, collinear_tol1e-9) print(fNumber of lines with 4 points: {num_rows}) print_points_and_lines(test_points, lines, (0,10))运行这个测试你可能看到0条或少数几条线因为随机点很难恰好4点共线。这正说明了问题的难度也说明我们需要优化算法来主动“制造”共线。4. 优化版本实现哈希法与模拟退火搜索基础版本只能评估随机布局的好坏。为了找到更好的布局我们需要一个搜索算法。同时评估函数也需要优化以支持频繁调用。我们先优化评估函数再实现模拟退火。4.1 高效评估基于直线哈希的检测我们实现之前提到的哈希法。思路是枚举所有点对计算它们确定的直线将直线标准化为一个可哈希的键并记录该直线上的所有点。def line_representation(p1: Point, p2: Point, tolerance: float 1e-10): 给定两点返回一个代表其所确定直线的标准化元组可用于哈希。 处理垂直线斜率无穷大和浮点数精度。 方法使用直线的一般式 Ax By C 0并标准化 (A, B, C)。 标准化规则使得 A, B, C 的绝对值尽可能小且 C 0如果 C0 则 B 0如果 B0 则 A 0。 x1, y1 p1 x2, y2 p2 # 计算向量差 dx x2 - x1 dy y2 - y1 # 如果两点非常接近视为同一点返回None忽略 if abs(dx) tolerance and abs(dy) tolerance: return None # 直线一般式参数A dy, B -dx, C dx*y1 - dy*x1 A dy B -dx C dx * y1 - dy * x1 # 标准化除以最大公约数浮点数近似并确保第一个非零系数为正 # 由于是浮点数我们用一个比例因子归一化使A和B的平方和为1单位法向量 norm math.hypot(A, B) if norm tolerance: # 这种情况理论上不会发生除非dx和dy都为零已处理 return None A, B, C A / norm, B / norm, C / norm # 确保“标准形式”如果 C 0或者 C 0 且 B 0或者 C 0 且 B 0 且 A 0 # 否则全体取反 if C -tolerance or (abs(C) tolerance and B -tolerance) or (abs(C) tolerance and abs(B) tolerance and A -tolerance): A, B, C -A, -B, -C # 为了哈希将参数四舍五入到一定精度以处理浮点误差 precision 12 # 小数点后12位 key (round(A, precision), round(B, precision), round(C, precision)) return key def evaluate_hash(points: List[Point], collinear_tol: float 1e-10) - Tuple[int, List[Set[int]]]: 使用哈希法评估找出所有包含至少4个点的直线。 Args: points: 点列表。 collinear_tol: 用于判断两点是否过近的容差。 Returns: (行数, 直线列表) 每条直线是一个点索引集合。 n len(points) line_map {} # key - set of point indices for i in range(n): for j in range(i1, n): key line_representation(points[i], points[j], collinear_tol) if key is None: continue # 忽略重合或过近的点对 if key not in line_map: line_map[key] set() line_map[key].add(i) line_map[key].add(j) # 筛选出点数4的直线 valid_lines [indices for indices in line_map.values() if len(indices) 4] num_rows len(valid_lines) return num_rows, valid_lines这个evaluate_hash函数比暴力枚举快得多O(n^2) vs O(n^4)并且能自然地处理一条直线上有任意多个点的情况。它直接给出了每条直线上的所有点索引。4.2 模拟退火算法实现现在我们实现模拟退火算法来搜索最优布局。def simulated_annealing(initial_points: List[Point], iterations: int 50000, start_temp: float 10.0, end_temp: float 0.01, step_size: float 2.0, eval_func evaluate_hash) - Tuple[List[Point], int, List[Set[int]]]: 模拟退火算法优化点布局以最大化共线直线数。 Args: initial_points: 初始点集。 iterations: 总迭代次数。 start_temp: 初始温度。 end_temp: 终止温度。 step_size: 点扰动的最大步长在坐标每个维度上。 eval_func: 评估函数默认为evaluate_hash。 Returns: (best_points, best_score, best_lines): 最优的点集、得分和直线信息。 current_points initial_points[:] # 深拷贝 current_score, _ eval_func(current_points) best_points current_points[:] best_score current_score best_lines [] n len(current_points) temp start_temp decay (end_temp / start_temp) ** (1.0 / iterations) for iter in range(iterations): # 1. 生成邻居随机扰动一个点 new_points current_points[:] idx random.randrange(n) x, y new_points[idx] # 在当前点附近随机扰动 new_x x random.uniform(-step_size, step_size) new_y y random.uniform(-step_size, step_size) # 可选将坐标限制在边界内如[0,100] new_x max(0, min(100, new_x)) new_y max(0, min(100, new_y)) new_points[idx] (new_x, new_y) # 2. 评估新解 new_score, new_lines eval_func(new_points) # 3. 决定是否接受新解 delta new_score - current_score if delta 0: # 新解更好总是接受 accept True else: # 新解更差以一定概率接受 prob math.exp(delta / temp) # 注意delta为负 accept random.random() prob if accept: current_points new_points current_score new_score if current_score best_score: best_points current_points[:] best_score current_score best_lines new_lines print(fIteration {iter}: New best score {best_score}, temp{temp:.4f}) # 4. 降温 temp * decay return best_points, best_score, best_lines4.3 随机重启框架单次模拟退火很可能陷入局部最优。我们需要多次随机重启保留历史最佳。def random_restart_search(num_restarts: int 50, points_per_restart: int 20, sa_iterations: int 20000, **sa_kwargs) - Tuple[List[Point], int, List[Set[int]]]: 随机重启搜索框架。 Args: num_restarts: 重启次数。 points_per_restart: 每次重启生成的点数。 sa_iterations: 每次模拟退火的迭代次数。 **sa_kwargs: 传递给simulated_annealing的其他参数。 Returns: 全局最优的点集、得分和直线信息。 global_best_points None global_best_score -1 global_best_lines [] for restart in range(num_restarts): print(f\n--- Random Restart {restart1}/{num_restarts} ---) init_points generate_random_points(points_per_restart, (0, 100)) best_points, best_score, best_lines simulated_annealing( init_points, iterationssa_iterations, **sa_kwargs ) if best_score global_best_score: global_best_score best_score global_best_points best_points[:] global_best_lines best_lines print(f - New global best score: {global_best_score}) print(f\n Global Best Score: {global_best_score} ) return global_best_points, global_best_score, global_best_lines4.4 运行优化与结果分析现在让我们运行这个优化程序看看能找到多少条线。if __name__ __main__: # 为了节省时间我们先运行一个较小规模的搜索 final_points, final_score, final_lines random_restart_search( num_restarts20, points_per_restart20, sa_iterations10000, start_temp5.0, end_temp0.1, step_size5.0, eval_funcevaluate_hash # 使用高效的哈希评估 ) print(f\nFinal layout (first 10 points):) for i, pt in enumerate(final_points[:10]): print(f {i}: ({pt[0]:.2f}, {pt[1]:.2f})) print(f\nLines found ({final_score} lines):) for i, line_indices in enumerate(final_lines): points_on_line [final_points[idx] for idx in line_indices] # 简单计算一下直线的斜率用于展示 if len(points_on_line) 2: p1, p2 points_on_line[0], points_on_line[1] dx p2[0] - p1[0] dy p2[1] - p1[1] if abs(dx) 1e-10: slope dy / dx print(f Line {i}: {len(line_indices)} points, slope ~ {slope:.2f}, points indices: {sorted(line_indices)}) else: print(f Line {i}: {len(line_indices)} points, vertical line, points indices: {sorted(line_indices)})运行这个程序可能需要一些时间取决于迭代次数和重启次数。在我的测试中经过数万次迭代和多次重启程序通常能找到包含5到8条满足条件的直线的布局。这已经比完全随机布局通常0-1条好得多。实操心得模拟退火的参数对结果影响很大。start_temp初始温度不宜过低否则算法退化成爬山法容易陷入局部最优。step_size扰动步长也很关键步长太大搜索过于随机步长太小难以跳出当前区域。我通常从step_size5.0相对于坐标范围0-100开始并在搜索后期逐步减小。iterations迭代次数需要足够多让温度能缓慢降低到end_temp。一个常见的策略是先进行几轮“粗搜索”参数范围大迭代次数少找到有希望的区域再进行“精搜索”缩小步长增加迭代。5. 高级优化与问题深潜通过基本的模拟退火我们可能找到了一个包含若干条4点共线直线的布局。但“20棵树每行4棵”的已知理论最大值是多少呢这是一个著名的数学问题。实际上对于n棵树每行k棵最大行数问题被称为“果园问题”。对于20棵树每行4棵已知的最优解是18条直线。这是一个非常精巧的几何构造与我们随机搜索得到的结果相去甚远。这说明我们的简单随机扰动策略很难发现这种高度对称和结构化的最优解。5.1 为何简单模拟退火难以找到最优解搜索空间巨大且非凸20个点在平面上的布局空间是40维的连续空间。最优解18条线位于这个空间中的一个非常狭窄的“尖峰”上。随机游走很难恰好命中这个区域。目标函数不平滑移动一个点可能会同时创建或破坏多条直线。目标函数直线数量是离散的且变化不连续这给基于梯度的优化模拟退火是随机梯度带来困难。需要高度对称性最优解往往具有特定的对称性如正五边形与星形结构的组合或者基于射影几何的配置。单纯的局部扰动无法产生这种全局对称结构。5.2 改进策略探索为了向最优解靠近我们需要更智能的搜索策略启发式初始化不要完全随机初始化。可以尝试从一些已知的、具有高共线潜力的图案开始例如将点放在一个规则网格上但网格本身可能不会产生4点斜线。将点放在几个交于一点的直线上增加共线机会。使用已知的较小规模最优解如9棵树每行3棵的“三三排列”作为子结构。针对性扰动吸引力对于已经接近共线的三个点可以轻微移动第四个点使其精确共线。排斥力对于密集区域可以施加轻微排斥避免点堆积让布局更均匀可能发现新的共线机会。直线对齐操作随机选择一条已有3个点的直线尝试移动其他点使其落到这条直线上。混合算法将模拟退火与局部搜索结合。在模拟退火接受一个新解后可以对其进行一轮快速的局部贪心优化例如尝试微调每个点看能否立即增加直线数。使用更专业的优化库对于这类连续优化问题可以尝试使用像SciPy中的优化算法如差分进化、盆地跳跃等它们可能比我们自写的模拟退火更强大。引入几何约束如果我们知道最优解可能与某些几何图形如五角星、正多边形有关可以在搜索中引入偏置让点更倾向于落在这些图形的顶点或交点上。5.3 实现一个简单的“对齐”扰动算子让我们在模拟退火中增加一种新的邻居生成方式以一定概率执行“对齐”操作而不是纯粹随机扰动。def neighbor_with_alignment(points: List[Point], step_size: float, align_prob: float 0.3): 生成邻居点集。以align_prob的概率尝试执行“对齐”操作否则进行普通随机扰动。 对齐操作随机选择一条已有至少3个点的直线再随机选择一个不在这条直线上的点 将该点移动到这条直线上在直线上随机选择一个位置。 n len(points) new_points points[:] if random.random() align_prob: # 尝试对齐操作 _, lines evaluate_hash(points) # 获取当前所有直线点数3? 这里evaluate_hash只返回4的我们需要3的 # 我们需要一个函数来获取所有至少3点共线的直线 # 为了简化我们临时计算所有3点共线的情况或者复用line_map但筛选点数3的直线 # 这里我们简单地从lines中选点数4如果lines为空则回退到随机扰动 if lines: line random.choice(lines) if len(line) 3: # 选择这条直线上的两个点来确定直线 line_indices list(line) p1_idx, p2_idx random.sample(line_indices, 2) p1, p2 points[p1_idx], points[p2_idx] # 选择一个不在这条直线上的点 all_indices set(range(n)) off_line_indices list(all_indices - line) if off_line_indices: move_idx random.choice(off_line_indices) # 计算直线上一个随机位置 # 直线参数方程: P p1 t * (p2 - p1) t random.uniform(-2.0, 3.0) # 可以扩展到线段外 new_x p1[0] t * (p2[0] - p1[0]) new_y p1[1] t * (p2[1] - p1[1]) # 限制边界 new_x max(0, min(100, new_x)) new_y max(0, min(100, new_y)) new_points[move_idx] (new_x, new_y) return new_points # 普通随机扰动 idx random.randrange(n) x, y new_points[idx] new_x x random.uniform(-step_size, step_size) new_y y random.uniform(-step_size, step_size) new_x max(0, min(100, new_x)) new_y max(0, min(100, new_y)) new_points[idx] (new_x, new_y) return new_points然后修改simulated_annealing函数用neighbor_with_alignment替换原来的邻居生成部分。这个操作增加了搜索的导向性可能会更快地增加直线数量。6. 结果验证、可视化与扩展思考经过长时间运行例如数百次重启每次数万次迭代我们的程序可能可以稳定地找到包含10条左右直线的布局。虽然离理论最大值18还有距离但这已经是一个显著的进步并且完美地展示了计算模拟在解决复杂组合优化问题上的能力。6.1 结果验证与输出我们可以将找到的最佳布局保存下来并用更直观的方式可视化。def save_layout(points: List[Point], filename: str best_layout.txt): 将点坐标保存到文件 with open(filename, w) as f: for i, (x, y) in enumerate(points): f.write(f{i}: {x:.6f}, {y:.6f}\n) print(fLayout saved to {filename}) def visualize_with_matplotlib(points: List[Point], lines: List[Set[int]]): 使用matplotlib可视化点和直线。 需要安装matplotlib: pip install matplotlib try: import matplotlib.pyplot as plt except ImportError: print(Matplotlib not installed. Install with pip install matplotlib to visualize.) return plt.figure(figsize(10, 10)) # 绘制点 xs, ys zip(*points) plt.scatter(xs, ys, colorblue, s100, zorder5) for i, (x, y) in enumerate(points): plt.text(x, y, str(i), fontsize12, hacenter, vacenter, colorwhite) # 绘制直线 colors plt.cm.tab10.colors for i, line_indices in enumerate(lines): if len(line_indices) 2: continue # 获取这条直线上的所有点 line_points [points[idx] for idx in line_indices] # 为了画线我们需要找到这条直线的两个端点最远的两个点 # 简单起见我们计算点的凸包或者直接取最小和最大的x对应的点如果直线不垂直 lxs, lys zip(*line_points) # 拟合一条直线最小二乘来绘制 import numpy as np x_arr np.array(lxs) y_arr np.array(lys) A np.vstack([x_arr, np.ones(len(x_arr))]).T m, c np.linalg.lstsq(A, y_arr, rcondNone)[0] # 生成直线的x范围 x_min, x_max min(x_arr), max(x_arr) x_plot np.array([x_min - 5, x_max 5]) # 延长一些 y_plot m * x_plot c color colors[i % len(colors)] plt.plot(x_plot, y_plot, colorcolor, linewidth2, alpha0.7, labelfLine {i} ({len(line_indices)} pts)) plt.xlim(-5, 105) plt.ylim(-5, 105) plt.gca().set_aspect(equal, adjustablebox) plt.title(fBest Layout with {len(lines)} Lines) plt.legend(locbest) plt.grid(True, alpha0.3) plt.show() # 在主程序找到最优解后调用 if __name__ __main__: # ... 运行搜索算法获取 final_points, final_score, final_lines ... save_layout(final_points, best_20trees.txt) visualize_with_matplotlib(final_points, final_lines)6.2 已知最优解与数学背景正如前文所述20棵树每行4棵的最大行数问题是著名的果园问题的一个实例。已知的经典解是通过构造一个5阶射影平面上的配置或者利用五角星和正五边形的重叠结构来实现的。这种结构具有高度的对称性20个点树位于5条同心五角星和正五边形的交点以及中心点上从而产生了18条直线每条直线上恰好有4个点。我们的模拟退火算法是一种元启发式算法它擅长在复杂的搜索空间中寻找较好的解但无法保证找到数学上的最优解。要找到18条直线的精确解通常需要利用问题的组合几何性质进行构造性证明或者使用更专门的符号计算和约束求解工具。6.3 项目扩展方向这个项目虽然围绕一个具体问题展开但其方法可以扩展到许多其他领域推广问题可以很容易地修改代码研究“n棵树每行k棵”的最大行数问题例如经典的“9棵树种10行每行3棵”。算法比较可以在这个问题上对比不同优化算法的性能如遗传算法、粒子群优化、蚁群算法等看看哪种元启发式算法更适合此类几何排列问题。约束优化引入更多约束例如树与树之间最小距离模拟实际种植或者点必须落在整数网格上离散优化。交互式探索构建一个图形界面允许用户手动拖动点实时计算并显示共线情况结合算法建议进行人机协同探索。连接理论数学将找到的近似最优解与已知的射影几何配置进行对比尝试理解其结构甚至可以尝试用算法来“发现”已知的最优图案。通过这个项目我们不仅练习了Python编程、几何计算、算法设计模拟退火和问题建模更重要的是体验了如何用计算思维去探索和逼近一个数学难题。它完美地展示了编程作为“实验数学”工具的强大之处。尽管我们可能无法用简单的模拟退火敲开18条直线的大门但在这个过程中学到的迭代优化、问题分解和调试技巧无疑是宝贵的财富。下次当你遇到一个看似无从下手的排列组合或几何极值问题时不妨想想能不能写个程序让计算机帮我找找答案