I was inspired to write a solver for this kind of puzzle (it's not pretty but it seems to work well): import numpy as np import z3 import matplotli

I was inspired to write a solver for this kind of puzzle (it's not pretty but it seems to work well): import numpy as np import z3 import matplotlib.pyplot as plt from sklearn.svm import LinearSVC from sklearn.cluster import KMeans def normalize(v): return v / np.sqrt(np.dot(v, v)) def line_normal(v): return np.array([v[1], -v[0]]) def line_segment_intersection_test(e0, e1): a, c, b, d = e0[0], e1[0], e0[1] - e0[0], e1[1] - e1[0] n = line_normal(b) l = np.dot(n, d) if l == 0: return False t1 = np.dot(n, a - c) / l if not 0 <= t1 <= 1: return False if b[0] == 0: return False t0 = (c[0] + t1 * d[0] - a[0]) / b[0] return 0 <= t0 <= 1 def min_cover(n, subsets): #some simple pruning makes z3 run much faster subsets = list(enumerate(map(set, subsets))) subsets.sort(key=lambda t: len(t[1])) redundant = set([]) for i, (_, s) in enumerate(subsets): for j in range(i + 1, len(subsets)): _, s1 = subsets[j] if all(x in s1 for x in s): redundant.add(i) break subsets = [t for i, t in enumerate(subsets) if i not in redundant] s = z3.Solver() included_subsets = [z3.Bool('subset_%d_included' % i) for i, _ in subsets] for i in range(n): covering_subsets = [j for j, (_, s) in enumerate(subsets) if i in s] s.add(z3.Or(*[included_subsets[j] for j in covering_subsets])) for bound in range(len(subsets)+1): print("trying %d cut%s..." % (bound, 's' * (bound != 1))) s.push() s.add(z3.PbEq([(v, 1) for v in included_subsets], bound)) if s.check() == z3.sat: model = s.model() return [subsets[i][0] for i, v in enumerate(included_subsets) if model.eval(v)] s.pop() return [] def solve(islands): edges = [] islands = [[np.array(v) for v in island] for island in islands] print("enumerating edges...") for island in islands: for i, v_i in enumerate(island): for j, v_j in enumerate(island): if j == i: break edges.append((v_i, v_j)) vertices = np.array([v for island in islands for v in island]) x, y = vertices.T x_min, y_min = np.min(vertices, axis=0) x_mean, y_mean = np.mean(vertices, axis=0) x_max, y_max = np.max(vertices, axis=0) d_max = 2 * np.sqrt((y_max - y_min) ** 2 + (x_max - x_min) ** 2) cuts = [] print("enumerating cuts...") for i, v_i in enumerate(vertices): for j, v_j in enumerate(vertices): if j == i: break dv = normalize(v_j - v_i) n = line_normal(dv) eps = 1e-7 #try both cut orientations (left+right and right+left), and extend the cuts to poke out of the domain c0 = (v_i + eps * n - d_max * dv, v_j - eps * n + d_max * dv) c1 = (v_i - eps * n - d_max * dv, v_j + eps * n + d_max * dv) cuts.append(c0) cuts.append(c1) print("finding intersections...") cut_sets = [[] for _ in cuts] for i, cut in enumerate(cuts): for j, edge in enumerate(edges): if line_segment_intersection_test(cut, edge): cut_sets[i].append(j) print("solving for min-cover...") included_cuts = [cuts[i] for i in min_cover(len(edges), cut_sets)] print("refining cuts using SVM...") refined_cuts = [] for cut in included_cuts: v0, v1 = cut dv = v1 - v0 n = line_normal(dv) side = (np.sum((vertices - v0) * n, axis=-1) > 0).astype(np.long) svm = LinearSVC(C=1e9, loss='hinge', max_iter=10000) svm.fit(vertices, side) if svm.score(vertices, side) == 1: m, c = svm.coef_[0], svm.intercept_[0] dv = line_normal(m) if m[1] == 0: v0 = np.array([-c/m[0], y_mean]) else: v0 = np.array([x_mean, -(c+m[0] * x_mean)/m[1]]) else: print("SVM fitting failed: (falling back to grazing cut)") dv = normalize(dv) proj = np.sum((vertices - v0) * dv, axis=-1) proj_min, proj_max = np.min(proj), np.max(proj) proj_range = proj_max - proj_min refined_cuts.append([v0 + (proj_min - 0.1 * proj_range) * dv, v0 + (proj_max + 0.1 * proj_range) * dv]) return refined_cuts vertices = np.random.uniform(-1, 1, (15, 2)) n_clusters = 5 k_means = KMeans(n_clusters=n_clusters).fit(vertices) labels = k_means.labels_ islands = [[v for j, v in enumerate(vertices) if labels[j] == l] for l in range(n_clusters)] cuts = solve(islands) for island in islands: plt.scatter(*np.array(island).T, marker='+') for cut in cuts: plt.plot(*np.array(cut).T, color='black') plt.show()