Skip to content

Improve performance and structure of polygon packing algorithm #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Samt4y wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Improve performance and structure of polygon packing algorithm #1

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
292 changes: 133 additions & 159 deletions polygon_packer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,195 +1,169 @@
import numpy as np
from scipy.optimize import basinhopping, minimize
import matplotlib.pyplot as ppt
import matplotlib.pyplot as plt
from numba import njit
from joblib import Parallel, delayed
import argparse
from io import BytesIO
import requests

arg_parser = argparse.ArgumentParser()
arg_parser.add_argument("inner_polygons", type=int, help="Number of inner polygons")
arg_parser.add_argument("inner_sides", type=int, help="Number of sides of the inner polygons")
arg_parser.add_argument("container_sides", type=int, help="Number of sides of the container polygon")
arg_parser.add_argument("--attempts", type=int, default=1000, help="Number of attempts to run")
arg_parser.add_argument("--tolerance", type=float, default=1e-8, help="Overlap penalty tolerance. Probably best left at default")
arg_parser.add_argument("--finalstep", type=float, default=0.0001, help="How small the last theoretical step in container size decrease will be (it gets smaller over time)")
args = arg_parser.parse_args()

parser = argparse.ArgumentParser()
parser.add_argument("inner_polygons", type=int)
parser.add_argument("inner_sides", type=int)
parser.add_argument("container_sides", type=int)
parser.add_argument("--attempts", type=int, default=1000)
parser.add_argument("--tolerance", type=float, default=1e-8)
parser.add_argument("--finalstep", type=float, default=1e-4)
args = parser.parse_args()

N = args.inner_polygons
nsi = args.inner_sides
nsc = args.container_sides
attempts = args.attempts
penalty_tolerance = args.tolerance
final_step_size = args.finalstep
penalty_tol = args.tolerance
final_step = args.finalstep

def regular_polygon(n):
angles = np.linspace(0, 2*np.pi, n, endpoint=False)
verts = np.column_stack((np.cos(angles), np.sin(angles)))
normals = np.column_stack((np.cos(angles + np.pi/n), np.sin(angles + np.pi/n)))
return verts, normals

unit_polygon_angles = np.linspace(0, 2 * np.pi, nsi, endpoint=False)
unit_polygon_vertices = np.column_stack((np.cos(unit_polygon_angles), np.sin(unit_polygon_angles)))
unit_polygon_vectors = np.column_stack((np.cos(unit_polygon_angles + np.pi / nsi), np.sin(unit_polygon_angles + np.pi / nsi)))
unit_container_angles = np.linspace(0, 2 * np.pi, nsc, endpoint=False)
unit_container_vertices = np.column_stack((np.cos(unit_container_angles), np.sin(unit_container_angles)))
unit_container_vectors = np.column_stack((np.cos(unit_container_angles + np.pi / nsc), np.sin(unit_container_angles + np.pi / nsc)))
unit_container_apothem = np.cos(np.pi / nsc)
unit_poly, unit_poly_normals = regular_polygon(nsi)
unit_cont, unit_cont_normals = regular_polygon(nsc)
container_apothem = np.cos(np.pi / nsc)

@njit(cache=True)
def transform_polygon(x, y, a, vertices):
n_vertices = vertices.shape[0]
transformed = np.empty_like(vertices)
for i in range(n_vertices):
vx = vertices[i, 0]
vy = vertices[i, 1]
transformed[i, 0] = x + (vx * np.cos(a) - vy * np.sin(a))
transformed[i, 1] = y + (vx * np.sin(a) + vy * np.cos(a))
return transformed
def transform(x, y, angle, verts):
c, s = np.cos(angle), np.sin(angle)
out = np.empty_like(verts)
for i in range(verts.shape[0]):
vx, vy = verts[i]
out[i, 0] = x + vx*c - vy*s
out[i, 1] = y + vx*s + vy*c
return out

@njit(cache=True)
def rotate_vectors(a, vectors):
n_vectors = vectors.shape[0]
rotated = np.empty_like(vectors)
for i in range(n_vectors):
vecx = vectors[i, 0]
vecy = vectors[i, 1]
rotated[i, 0] = vecx * np.cos(a) - vecy * np.sin(a)
rotated[i, 1] = vecx * np.sin(a) + vecy * np.cos(a)
return rotated

def rotate(angle, vecs):
c, s = np.cos(angle), np.sin(angle)
out = np.empty_like(vecs)
for i in range(vecs.shape[0]):
vx, vy = vecs[i]
out[i, 0] = vx*c - vy*s
out[i, 1] = vx*s + vy*c
return out

@njit(cache=True)
def poking_penalty(vertices, S):
def container_penalty(verts, S):
limit = container_apothem * S
penalty = 0.0
limit = unit_container_apothem * S
for v in range(vertices.shape[0]):
vx = vertices[v, 0]
vy = vertices[v, 1]
for v in range(verts.shape[0]):
x, y = verts[v]
for i in range(nsc):
distance = vx * unit_container_vectors[i, 0] + vy * unit_container_vectors[i, 1]
if distance > limit:
diff = distance - limit
d = x * unit_cont_normals[i, 0] + y * unit_cont_normals[i, 1]
if d > limit:
diff = d - limit
penalty += diff * diff
return penalty

@njit(cache=True)
def bh_function(values, S):
def sat_overlap(poly1, poly2, axes):
min_overlap = 1e18
for ax in axes:
axx, axy = ax

min1 = 1e18
max1 = -1e18
for v in poly1:
d = v[0]*axx + v[1]*axy
if d < min1: min1 = d
if d > max1: max1 = d

min2 = 1e18
max2 = -1e18
for v in poly2:
d = v[0]*axx + v[1]*axy
if d < min2: min2 = d
if d > max2: max2 = d

overlap = min(max1, max2) - max(min1, min2)
if overlap <= 0:
return 0.0 # no collision
if overlap < min_overlap:
min_overlap = overlap

return min_overlap

@njit(cache=True)
def objective(vals, S):
penalty = 0.0
polygon_array = np.zeros((N, nsi, 2))
vector_array = np.zeros((N, nsi, 2))
for i in range(N):
posx = values[i * 3]
posy = values[i * 3 + 1]
rot = values[i * 3 + 2]
polygon_array[i] = transform_polygon(posx, posy, rot, unit_polygon_vertices)
vector_array[i] = rotate_vectors(rot, unit_polygon_vectors)

penalty += poking_penalty(polygon_array[i], S)
polys = np.empty((N, nsi, 2))
axes = np.empty((N, nsi, 2))

for i in range(N):
for j in range(i + 1, N):
collision = True
min_overlap = 100000000000000000000.0
for vec in range(nsi * 2):
if vec < nsi:
x_axis = vector_array[i][vec, 0]
y_axis = vector_array[i][vec, 1]
else:
x_axis = vector_array[j][vec - nsi, 0]
y_axis = vector_array[j][vec - nsi, 1]

min_1 = 100000000000000000000.0
max_1 = -100000000000000000000.0
for vert in range(nsi):
dotp = polygon_array[i][vert, 0] * x_axis + polygon_array[i][vert, 1] * y_axis
if dotp < min_1:
min_1 = dotp
if dotp > max_1:
max_1 = dotp

min_2 = 100000000000000000000.0
max_2 = -100000000000000000000.0
for vert in range(nsi):
dotp = polygon_array[j][vert, 0] * x_axis + polygon_array[j][vert, 1] * y_axis
if dotp < min_2:
min_2 = dotp
if dotp > max_2:
max_2 = dotp

overlap = min(max_1, max_2) - max(min_1, min_2)
if overlap <= 0:
collision = False
break
if overlap < min_overlap:
min_overlap = overlap

if collision:
penalty += min_overlap * min_overlap

return penalty
x, y, a = vals[i*3:i*3+3]
polys[i] = transform(x, y, a, unit_poly)
axes[i] = rotate(a, unit_poly_normals)
penalty += container_penalty(polys[i], S)

def repetition(seed):
print("Attempt", seed)
for i in range(N):
for j in range(i+1, N):
combined_axes = np.vstack((axes[i], axes[j]))
overlap = sat_overlap(polys[i], polys[j], combined_axes)
if overlap > 0:
penalty += overlap * overlap

return penalty

def run_attempt(seed):
np.random.seed(seed)
dynamic_S = np.sqrt(N) * (2 + np.random.rand() * 2)
initial_S = dynamic_S

if np.random.rand() < 0.5:
x0 = np.random.uniform(-dynamic_S/2, dynamic_S/2, N * 3)
else:
grid_linspace = np.linspace(-dynamic_S/2 * 0.9, dynamic_S/2 * 0.9, int(np.ceil(np.sqrt(N))))
xx, yy = np.meshgrid(grid_linspace, grid_linspace)
grid_points = np.column_stack((xx.flatten(), yy.flatten()))[:N]
x0 = np.zeros(N * 3)
x0[0::3] = grid_points[:, 0]
x0[1::3] = grid_points[:, 1]
x0[2::3] = np.random.uniform(0, 2 * np.pi, N)

last_valid_x = x0.copy()
last_valid_S = dynamic_S

S = np.sqrt(N) * (2 + np.random.rand()*2)
initial_S = S

x0 = np.random.uniform(-S/2, S/2, N*3)

best_x = x0.copy()
best_S = S

while True:
minimized = minimize(bh_function, x0, args=(dynamic_S,), method="L-BFGS-B", tol=1e-8)
multiplier = 0.9999 - (dynamic_S - np.sqrt(N) * nsi / nsc) * 0.0099 / (initial_S - np.sqrt(N) * nsi / nsc)
multiplier = 1 - final_step_size - (dynamic_S - np.sqrt(N) * nsi / nsc) * (0.01 - final_step_size) / (initial_S - np.sqrt(N) * nsi / nsc)
if minimized.fun < penalty_tolerance:
last_valid_x = minimized.x.copy()
last_valid_S = dynamic_S.copy()
x0 = minimized.x * multiplier
dynamic_S *= multiplier
res = minimize(objective, x0, args=(S,), method="L-BFGS-B")

shrink = 1 - final_step - (S - np.sqrt(N)*nsi/nsc) * (0.01 - final_step) / (initial_S - np.sqrt(N)*nsi/nsc)

if res.fun < penalty_tol:
best_x, best_S = res.x.copy(), S
x0 = res.x * shrink
S *= shrink
else:
bh_result = basinhopping(
lambda x, s: bh_function(x, s),
x0,
minimizer_kwargs={'method': 'L-BFGS-B', 'args': (dynamic_S,), 'tol': 1e-8},
niter=50,
T=0.1,
stepsize=0.1
)
if bh_result.fun < penalty_tolerance:
last_valid_x = bh_result.x.copy()
last_valid_S = dynamic_S
x0 = bh_result.x * multiplier
dynamic_S *= multiplier
bh = basinhopping(objective, x0,
minimizer_kwargs={"args": (S,)},
niter=50)
if bh.fun < penalty_tol:
best_x, best_S = bh.x.copy(), S
x0 = bh.x * shrink
S *= shrink
else:
break
return last_valid_S, last_valid_x

best_S = float("inf")
best_values = None
results = Parallel(n_jobs=-1, prefer="processes")(delayed(repetition)(i) for i in range(attempts))

for s, values in results:
if s < best_S:
best_S = s
best_values = values

print("Final side length:", best_S * np.sin(np.pi / nsc) / np.sin(np.pi / nsi))

final_positions = positions = best_values.reshape((N, 3))
fig, ax = ppt.subplots()
container_plot = np.vstack((unit_container_vertices * best_S, unit_container_vertices[0] * best_S))
ax.plot(container_plot[:,0], container_plot[:,1], color="#000000", linewidth=0.5)

return best_S, best_x

results = Parallel(n_jobs=-1)(delayed(run_attempt)(i) for i in range(attempts))
best_S, best_vals = min(results, key=lambda x: x[0])

side_len = best_S * np.sin(np.pi/nsc) / np.sin(np.pi/nsi)
print("Final side length:", side_len)

fig, ax = plt.subplots()

cont = np.vstack((unit_cont * best_S, unit_cont[0]*best_S))
ax.plot(cont[:,0], cont[:,1], 'k-', lw=0.5)

for i in range(N):
polygon = transform_polygon(best_values[i * 3], best_values[i * 3 + 1], best_values[i * 3 + 2], unit_polygon_vertices)
polygon_plot = np.vstack((polygon, polygon[0]))
ax.fill(polygon_plot[:,0], polygon_plot[:,1], "#CCCCCC", edgecolor="black", linewidth=0.5)
x, y, a = best_vals[i*3:i*3+3]
poly = transform(x, y, a, unit_poly)
poly = np.vstack((poly, poly[0]))
ax.fill(poly[:,0], poly[:,1], "#ccc", edgecolor="black", lw=0.5)

ax.set_aspect("equal")
ppt.title(f"Side length: {best_S * np.sin(np.pi / nsc) / np.sin(np.pi / nsi)}")
ppt.savefig(f"{N}_{nsi}_in_{nsc}.png")
plt.title(f"Side length: {side_len}")
plt.savefig(f"{N}_{nsi}_in_{nsc}.png")