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()