diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..d02e6ac
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+.venv
+polygon-packer-py/__pycache__
+polygon-packer-py/venv
+polygon-packer-rs/target
+polygon-packer-rs/Cargo.lock
+**/*.png
diff --git a/README.md b/README.md
index 17fe194..e695e03 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
This program can quickly solve the 2D bin packing problem for any number of any polygons inside any other polygon! It was the tool used to find all the optimal packings under the name "Ignacio Vallejo" on [Erich's Packing Center](https://erich-friedman.github.io/packing/).
-### How to use
+### How to use with Python
Download the python file, navigate to its location and run it like this:
`python3 polygon_packer.py [n] [nsi] [nsc]`
@@ -14,3 +14,24 @@ Optional parameters:
- `--attempts`: the total number of attempts to run. Increase to explore more possible packings. Defaults to 1000.
- `--tolerance`: the tolerance for the penalty function. More penalty reduces the margin of overlap but limits exporation. Defaults to an empirical sweetspot of 0.00000001.
- `--finalstep`: the container size is decreased by a smaller factor each time, to save compute at the beginning and achieve greater precision near the end. This sets the step size of the shrinkage which would correspond to the theoretical minimum container size (which, for most packings, will not actually be reached, so keep that in mind when setting this parameter). Defaults to 0.0001.
+
+### How to use with Rust
+
+- [Install The Rust Compiler.](https://rust-lang.org/learn/get-started/)
+
+- Clone this repository, change directory into `polygon-packer-rs`.
+
+- To compile with optimizations and run:
+
+```bash
+cargo run --release -- [n] [nsi] [nsc]
+```
+
+^ This may take a while.
+
+For further performance gains, set the environment variable `RUSTFLAGS` to
+`-C target-cpu=native`, and run the compiler again. Then, the compiler generates
+microarchitecture-specific machine code tailored for your CPU only.
+
+The part after `--release --` is the same arguments with the Python
+version.
\ No newline at end of file
diff --git a/polygon_packer.py b/polygon-packer-py/polygon_packer.py
similarity index 97%
rename from polygon_packer.py
rename to polygon-packer-py/polygon_packer.py
index 3fe2854..db94138 100644
--- a/polygon_packer.py
+++ b/polygon-packer-py/polygon_packer.py
@@ -1,193 +1,195 @@
-import numpy as np
-from scipy.optimize import basinhopping, minimize
-import matplotlib.pyplot as ppt
-from numba import njit
-from joblib import Parallel, delayed
-import argparse
-
-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()
-
-N = args.inner_polygons
-nsi = args.inner_sides
-nsc = args.container_sides
-attempts = args.attempts
-penalty_tolerance = args.tolerance
-final_step_size = args.finalstep
-
-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)
-
-@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
-
-@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
-
-@njit(cache=True)
-def poking_penalty(vertices, 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 i in range(nsc):
- distance = vx * unit_container_vectors[i, 0] + vy * unit_container_vectors[i, 1]
- if distance > limit:
- diff = distance - limit
- penalty += diff * diff
- return penalty
-
-@njit(cache=True)
-def bh_function(values, 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)
-
- 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
-
-def repetition(seed):
- print("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
-
- 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
- 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
- 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)
-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)
-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")
+import numpy as np
+from scipy.optimize import basinhopping, minimize
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as ppt
+from numba import njit
+from joblib import Parallel, delayed
+import argparse
+
+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()
+
+N = args.inner_polygons
+nsi = args.inner_sides
+nsc = args.container_sides
+attempts = args.attempts
+penalty_tolerance = args.tolerance
+final_step_size = args.finalstep
+
+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)
+
+@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
+
+@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
+
+@njit(cache=True)
+def poking_penalty(vertices, 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 i in range(nsc):
+ distance = vx * unit_container_vectors[i, 0] + vy * unit_container_vectors[i, 1]
+ if distance > limit:
+ diff = distance - limit
+ penalty += diff * diff
+ return penalty
+
+@njit(cache=True)
+def bh_function(values, 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)
+
+ 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
+
+def repetition(seed):
+ print("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
+
+ 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
+ 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
+ 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)
+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)
+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")
diff --git a/requirements.txt b/polygon-packer-py/requirements.txt
similarity index 100%
rename from requirements.txt
rename to polygon-packer-py/requirements.txt
diff --git a/polygon-packer-rs/Cargo.toml b/polygon-packer-rs/Cargo.toml
new file mode 100644
index 0000000..906465a
--- /dev/null
+++ b/polygon-packer-rs/Cargo.toml
@@ -0,0 +1,237 @@
+[package]
+name = "polygon-packer-rs"
+version = "0.2.0"
+edition = "2024"
+authors = ["Voxell Paladynee "]
+
+[dependencies]
+# bump allocator for temporary arena allocations to reduce transient heap use
+bumpalo = { version = "3", features = ["collections"] }
+# command line argument parsing
+clap = { version = "4.4", features = ["derive"] }
+# linear algebra library for geometric computations
+nalgebra = "0.34.2"
+# numeric library
+ndarray = "0.17.2"
+# optimization library
+nlopt = "0.6.0"
+# plot library
+plotters = "0.3"
+# parallelization
+rayon = "1.7"
+# Voxell's homegrown RNG library, used for pure-rust PCG family of random number
+# generators. https://www.pcg-random.org/index.html
+voxell_rng = "0.6.0"
+# Voxell's homegrown timer library
+voxell_timer = "1.2.2"
+
+[profile.release]
+# opt-level: 3 for full optimizations
+opt-level = 3
+# lto: allows library code to be more aggressively inlined
+# "thin" for a lighter mode, "fat" for max perf
+lto = "fat"
+# codegen-units: only 1 concurrent compiler so the optimizer
+# has maximum code context
+codegen-units = 1
+# panic: this potentially reduces instruction-cache pressure, but if it crashes
+# you won't get fancy errors.
+# set to "unwind" for fancy errors
+# set to "abort" for a very slight boost
+panic = "abort"
+
+# ============================================
+# Linter configuration. A linter is a program that pattern-matches on source
+# code to identify issues. I (@Paladynee) use a very strict linter
+# configuration.
+
+[lints.clippy]
+pedantic = { level = "warn", priority = -1 }
+cargo = { level = "allow", priority = -1 }
+restriction = { level = "warn", priority = -1 }
+nursery = { level = "warn", priority = -1 }
+
+mod_module_files = "deny"
+self_named_module_files = "allow"
+doc_paragraphs_missing_punctuation = "allow"
+allow_attributes = "allow"
+allow_attributes_without_reason = "allow"
+arithmetic_side_effects = "allow"
+as_conversions = "allow"
+as_pointer_underscore = "allow"
+as_underscore = "allow"
+assertions_on_result_states = "allow"
+big_endian_bytes = "allow"
+cognitive_complexity = "allow"
+arbitrary_source_item_ordering = "allow"
+dbg_macro = "allow"
+default_numeric_fallback = "allow"
+default_union_representation = "allow"
+deref_by_slicing = "allow"
+disallowed_script_idents = "allow"
+else_if_without_else = "allow"
+empty_drop = "allow"
+empty_enum_variants_with_brackets = "allow"
+empty_structs_with_brackets = "allow"
+exhaustive_enums = "allow"
+exhaustive_structs = "allow"
+expect_used = "allow"
+field_scoped_visibility_modifiers = "allow"
+float_arithmetic = "allow"
+fn_to_numeric_cast_any = "allow"
+host_endian_bytes = "allow"
+if_then_some_else_none = "allow"
+impl_trait_in_params = "allow"
+implicit_return = "allow"
+indexing_slicing = "allow"
+infinite_loop = "allow"
+inline_asm_x86_att_syntax = "allow"
+inline_asm_x86_intel_syntax = "allow"
+integer_division = "allow"
+integer_division_remainder_used = "allow"
+large_include_file = "allow"
+let_underscore_must_use = "allow"
+let_underscore_untyped = "allow"
+little_endian_bytes = "allow"
+lossy_float_literal = "allow"
+mem_forget = "allow"
+min_ident_chars = "allow"
+missing_assert_message = "allow"
+missing_asserts_for_indexing = "allow"
+missing_docs_in_private_items = "allow"
+missing_inline_in_public_items = "allow"
+missing_trait_methods = "allow"
+mixed_read_write_in_expression = "allow"
+module_name_repetitions = "allow"
+modulo_arithmetic = "allow"
+multiple_unsafe_ops_per_block = "allow"
+needless_raw_strings = "allow"
+non_ascii_literal = "allow"
+panic = "allow"
+panic_in_result_fn = "allow"
+partial_pub_fields = "allow"
+pointer_format = "allow"
+print_stderr = "allow"
+print_stdout = "allow"
+pub_use = "allow"
+pub_with_shorthand = "allow"
+pub_without_shorthand = "allow"
+question_mark_used = "allow"
+rc_buffer = "allow"
+redundant_type_annotations = "allow"
+ref_patterns = "allow"
+renamed_function_params = "allow"
+semicolon_inside_block = "allow"
+semicolon_outside_block = "allow"
+separated_literal_suffix = "allow"
+shadow_reuse = "allow"
+shadow_same = "allow"
+shadow_unrelated = "allow"
+single_call_fn = "allow"
+single_char_lifetime_names = "allow"
+std_instead_of_alloc = "allow"
+std_instead_of_core = "allow"
+str_to_string = "allow"
+string_add = "allow"
+string_slice = "allow"
+suspicious_xor_used_as_pow = "allow"
+tests_outside_test_module = "allow"
+todo = "allow"
+try_err = "allow"
+undocumented_unsafe_blocks = "allow"
+unimplemented = "allow"
+unnecessary_safety_comment = "allow"
+unnecessary_safety_doc = "allow"
+unnecessary_self_imports = "allow"
+unreachable = "allow"
+unseparated_literal_suffix = "allow"
+unused_result_ok = "allow"
+unused_trait_names = "allow"
+unwrap_in_result = "allow"
+unwrap_used = "allow"
+use_debug = "allow"
+wildcard_enum_match_arm = "allow"
+bool_to_int_with_if = "allow"
+borrow_as_ptr = "allow"
+cast_lossless = "allow"
+cast_possible_truncation = "allow"
+cast_precision_loss = "allow"
+cast_possible_wrap = "allow"
+cast_sign_loss = "allow"
+checked_conversions = "allow"
+default_trait_access = "allow"
+elidable_lifetime_names = "allow"
+empty_enums = "allow"
+enum_glob_use = "allow"
+explicit_deref_methods = "allow"
+explicit_into_iter_loop = "allow"
+explicit_iter_loop = "allow"
+if_not_else = "allow"
+ignore_without_reason = "allow"
+ignored_unit_patterns = "allow"
+index_refutable_slice = "allow"
+inline_always = "allow"
+into_iter_without_iter = "allow"
+items_after_statements = "allow"
+large_digit_groups = "allow"
+large_stack_arrays = "allow"
+large_types_passed_by_value = "allow"
+linkedlist = "allow"
+manual_assert = "allow"
+manual_is_power_of_two = "allow"
+many_single_char_names = "allow"
+match_bool = "allow"
+match_same_arms = "allow"
+match_wild_err_arm = "allow"
+match_wildcard_for_single_variants = "allow"
+missing_errors_doc = "allow"
+missing_fields_in_debug = "allow"
+missing_panics_doc = "allow"
+must_use_candidate = "allow"
+mut_mut = "allow"
+needless_bitwise_bool = "allow"
+needless_continue = "allow"
+needless_for_each = "allow"
+needless_pass_by_value = "allow"
+needless_raw_string_hashes = "allow"
+float_cmp = "allow"
+no_effect_underscore_binding = "allow"
+ptr_as_ptr = "allow"
+ptr_cast_constness = "allow"
+pub_underscore_fields = "allow"
+redundant_closure_for_method_calls = "allow"
+redundant_else = "allow"
+return_self_not_must_use = "allow"
+semicolon_if_nothing_returned = "allow"
+should_panic_without_expect = "allow"
+similar_names = "allow"
+single_match_else = "allow"
+string_add_assign = "allow"
+struct_excessive_bools = "allow"
+struct_field_names = "allow"
+too_many_lines = "allow"
+uninlined_format_args = "allow"
+unnecessary_box_returns = "allow"
+unnecessary_wraps = "allow"
+unnested_or_patterns = "allow"
+unreadable_literal = "allow"
+unused_self = "allow"
+used_underscore_binding = "allow"
+used_underscore_items = "allow"
+verbose_bit_mask = "allow"
+wildcard_imports = "allow"
+equatable_if_let = "allow"
+missing_const_for_fn = "allow"
+option_if_let_else = "allow"
+use_self = "allow"
+missing_safety_doc = "allow"
+manual_let_else = "allow"
+map_err_ignore = "allow"
+or_fun_call = "allow"
+blanket_clippy_restriction_lints = "allow"
+doc_markdown = "allow"
+decimal_literal_representation = "allow"
+
+[lints.rust]
+ambiguous_negative_literals = "warn"
+non_ascii_idents = "warn"
diff --git a/polygon-packer-rs/rustfmt.toml b/polygon-packer-rs/rustfmt.toml
new file mode 100644
index 0000000..95ea60c
--- /dev/null
+++ b/polygon-packer-rs/rustfmt.toml
@@ -0,0 +1,24 @@
+max_width = 80
+imports_granularity = "Item"
+edition = "2024"
+style_edition = "2024"
+group_imports = "StdExternalCrate"
+use_field_init_shorthand = true
+wrap_comments = true
+tab_spaces = 4
+newline_style = "Unix"
+reorder_imports = true
+reorder_impl_items = false
+reorder_modules = true
+format_code_in_doc_comments = true
+hard_tabs = false
+comment_width = 80
+format_macro_matchers = true
+imports_layout = "Vertical"
+merge_derives = true
+overflow_delimited_expr = true
+condense_wildcard_suffixes = true
+remove_nested_parens = true
+fn_params_layout = "Compressed"
+match_arm_blocks = false
+struct_lit_single_line = false
\ No newline at end of file
diff --git a/polygon-packer-rs/src/main.rs b/polygon-packer-rs/src/main.rs
new file mode 100644
index 0000000..2930820
--- /dev/null
+++ b/polygon-packer-rs/src/main.rs
@@ -0,0 +1,578 @@
+// changing `a * b + c` to `a.mul_add(b, c)` literally leads to more
+// verification errors, so we don't do that.
+#![allow(clippy::suboptimal_flops)]
+use std::f64::consts::PI;
+use std::iter;
+use std::time::Duration;
+
+use bumpalo::Bump;
+use clap::Parser;
+use nlopt::Algorithm;
+use nlopt::Nlopt;
+use nlopt::Target;
+use plotters::prelude::*;
+use rayon::prelude::*;
+use voxell_rng::rng::pcg_advanced::pcg_64::PcgInnerState64;
+use voxell_timer::time_fn;
+
+const MAX_ITERATIONS: usize = 5000;
+const BH_STEPS: usize = 50;
+const BH_STEP_SIZE: f64 = 0.2;
+const BH_TEMP: f64 = 0.1;
+const VERIFY_EPS: f64 = 1e-9;
+// const GRADIENT_EPS_REL: f64 = 1.0;
+const ABS_GRADIENT_STEP: f64 = 1e-8;
+// enable for C1-continuity issues
+const CENTRAL_DIFF: bool = true;
+
+#[derive(Parser, Debug, Clone)]
+pub struct Args {
+ /// Number of inner polygons
+ pub inner_polygons: usize,
+ /// Number of sides of the inner polygons
+ pub inner_sides: usize,
+ /// Number of sides of the container polygon
+ pub container_sides: usize,
+ /// Number of attempts to run
+ #[arg(long, default_value_t = 1000)]
+ pub attempts: usize,
+ /// Overlap penalty tolerance. Probably best left at default.
+ #[arg(long, default_value_t = 1e-8)]
+ pub tolerance: f64,
+ /// How small the last theoretical step in container size decrease will be
+ /// (it gets smaller over time)
+ #[arg(long, default_value_t = 0.0001)]
+ pub finalstep: f64,
+}
+
+fn main() {
+ eprintln!(
+ "Polygon Packer - Packing regular polygons into a regular container"
+ );
+ let args = Args::parse();
+ let geo = Geometry::new(args.inner_sides, args.container_sides);
+
+ println!(
+ "Starting packing with {} polygons ({} sides each) into a {}-sided container",
+ args.inner_polygons, args.inner_sides, args.container_sides
+ );
+ let (results, duration): (Vec<(f64, Vec)>, Duration) = time_fn(|| {
+ (0..args.attempts)
+ .into_par_iter()
+ .map(|i| repetition(i, &args, &geo))
+ .collect()
+ });
+
+ println!("Completed {} attempts in {:.2?}", args.attempts, duration);
+
+ if let Some((best_s, best_v)) = results
+ .into_iter()
+ .min_by(|a, b| a.0.partial_cmp(&b.0).unwrap())
+ {
+ let side_len = best_s * (PI / args.container_sides as f64).sin()
+ / (PI / args.inner_sides as f64).sin();
+ println!("Final side length: {}", side_len);
+
+ verify_results(&best_v, best_s, &args, &geo);
+ render(&args, &geo, best_s, &best_v);
+ }
+}
+
+struct Geometry {
+ inner_verts_x: Vec,
+ inner_verts_y: Vec,
+ inner_axes_x: Vec,
+ inner_axes_y: Vec,
+ cont_axes_x: Vec,
+ cont_axes_y: Vec,
+ apothem: f64,
+}
+
+impl Geometry {
+ fn new(nsi: usize, nsc: usize) -> Self {
+ let i_angle = 2.0 * PI / nsi as f64;
+ let c_angle = 2.0 * PI / nsc as f64;
+
+ let mut inner_verts_x = Vec::with_capacity(nsi);
+ let mut inner_verts_y = Vec::with_capacity(nsi);
+ for i in 0..nsi {
+ inner_verts_x.push((i as f64 * i_angle).cos());
+ inner_verts_y.push((i as f64 * i_angle).sin());
+ }
+
+ let mut inner_axes_x = Vec::with_capacity(nsi);
+ let mut inner_axes_y = Vec::with_capacity(nsi);
+ for i in 0..nsi {
+ let a = i as f64 * i_angle + PI / nsi as f64;
+ inner_axes_x.push(a.cos());
+ inner_axes_y.push(a.sin());
+ }
+
+ let mut cont_axes_x = Vec::with_capacity(nsc);
+ let mut cont_axes_y = Vec::with_capacity(nsc);
+ for i in 0..nsc {
+ let a = i as f64 * c_angle + PI / nsc as f64;
+ cont_axes_x.push(a.cos());
+ cont_axes_y.push(a.sin());
+ }
+
+ Self {
+ inner_verts_x,
+ inner_verts_y,
+ inner_axes_x,
+ inner_axes_y,
+ cont_axes_x,
+ cont_axes_y,
+ apothem: (PI / nsc as f64).cos(),
+ }
+ }
+}
+
+struct ScratchPad<'a> {
+ // (x, y) vertices for each polygon,
+ // flattened: [poly0_vx0, poly0_vy0, poly0_vx1, ...]
+ poly_verts: &'a mut [f64],
+ // (x, y) axes for each polygon, flattened
+ poly_axes: &'a mut [f64],
+}
+
+impl<'a> ScratchPad<'a> {
+ fn new(bump: &'a Bump, n_polys: usize, n_sides: usize) -> Self {
+ let v_size = n_polys * n_sides * 2;
+ let a_size = n_polys * n_sides * 2;
+ Self {
+ poly_verts: bump.alloc_slice_fill_copy(v_size, 0.0),
+ poly_axes: bump.alloc_slice_fill_copy(a_size, 0.0),
+ }
+ }
+}
+
+fn penalty(
+ v: &[f64], s: f64, args: &Args, geo: &Geometry, pad: &mut ScratchPad,
+) -> f64 {
+ let mut p = 0.0;
+ let limit = geo.apothem * s;
+ let nsi = args.inner_sides;
+
+ for i in 0..args.inner_polygons {
+ let (px, py, pr) = (v[i * 3], v[i * 3 + 1], v[i * 3 + 2]);
+ let (sin_r, cos_r) = pr.sin_cos();
+
+ for j in 0..nsi {
+ let vx = geo.inner_verts_x[j];
+ let vy = geo.inner_verts_y[j];
+ let tx = cos_r * vx - sin_r * vy + px;
+ let ty = sin_r * vx + cos_r * vy + py;
+
+ pad.poly_verts[(i * nsi + j) * 2] = tx;
+ pad.poly_verts[(i * nsi + j) * 2 + 1] = ty;
+
+ // for k in 0..geo.cont_axes_x.len() {
+ // let dist = tx * geo.cont_axes_x[k] + ty * geo.cont_axes_y[k];
+ // if dist > limit {
+ // p += (dist - limit).powi(2);
+ // }
+ // }
+ geo.cont_axes_x.iter().zip(&geo.cont_axes_y).for_each(
+ |(&ax, &ay)| {
+ let dist = tx * ax + ty * ay;
+ if dist > limit {
+ p += (dist - limit).powi(2);
+ }
+ },
+ );
+ }
+
+ for j in 0..nsi {
+ let ax = geo.inner_axes_x[j];
+ let ay = geo.inner_axes_y[j];
+ pad.poly_axes[(i * nsi + j) * 2] = cos_r * ax - sin_r * ay;
+ pad.poly_axes[(i * nsi + j) * 2 + 1] = sin_r * ax + cos_r * ay;
+ }
+ }
+
+ for i in 0..args.inner_polygons {
+ let i_off = i * nsi * 2;
+ for j in i + 1..args.inner_polygons {
+ let j_off = j * nsi * 2;
+ let mut min_overlap = 1e20;
+ let mut collision = true;
+
+ for k in 0..(nsi * 2) {
+ let (ax, ay) = if k < nsi {
+ (
+ pad.poly_axes[i_off + k * 2],
+ pad.poly_axes[i_off + k * 2 + 1],
+ )
+ } else {
+ let k2 = k - nsi;
+ (
+ pad.poly_axes[j_off + k2 * 2],
+ pad.poly_axes[j_off + k2 * 2 + 1],
+ )
+ };
+
+ let (min1, max1) =
+ project_soa(pad.poly_verts, i_off, nsi, ax, ay);
+ let (min2, max2) =
+ project_soa(pad.poly_verts, j_off, nsi, ax, ay);
+
+ let overlap = max1.min(max2) - min1.max(min2);
+ if overlap <= 0.0 {
+ collision = false;
+ break;
+ }
+ if overlap < min_overlap {
+ min_overlap = overlap;
+ }
+ }
+
+ if collision {
+ p += min_overlap * min_overlap;
+ }
+ }
+ }
+ p
+}
+
+#[inline(always)]
+fn project_soa(
+ verts: &[f64], offset: usize, n: usize, ax: f64, ay: f64,
+) -> (f64, f64) {
+ let mut min = f64::MAX;
+ let mut max = f64::MIN;
+ // for i in 0..n {
+ // let dot = verts[offset + i * 2] * ax + verts[offset + i * 2 + 1] *
+ // ay;
+ let relevant_verts = &verts[offset..offset + n * 2];
+ for i in 0..n {
+ let dot = relevant_verts[i * 2] * ax + relevant_verts[i * 2 + 1] * ay;
+ if dot < min {
+ min = dot;
+ }
+ if dot > max {
+ max = dot;
+ }
+ }
+ (min, max)
+}
+
+fn local_min(
+ x: &mut [f64], s: f64, args: &Args, geo: &Geometry, bump: &mut Bump,
+) -> Option> {
+ bump.reset();
+
+ struct OptCtx<'a> {
+ pad: ScratchPad<'a>,
+ x_tmp: Vec,
+ }
+
+ let ctx = OptCtx {
+ pad: ScratchPad::new(&*bump, args.inner_polygons, args.inner_sides),
+ x_tmp: x.to_vec(),
+ };
+
+ let obj = |xx: &[f64], grad: Option<&mut [f64]>, ctx: &mut OptCtx| -> f64 {
+ let f0 = penalty(xx, s, args, geo, &mut ctx.pad);
+ if let Some(g) = grad {
+ // let eps = f64::EPSILON.sqrt();
+ ctx.x_tmp.copy_from_slice(xx);
+ for i in 0..xx.len() {
+ // let step = eps * (GRADIENT_EPS_REL + xx[i].abs());
+ // ctx.x_tmp[i] += step;
+ // let f1 = penalty(&ctx.x_tmp, s, args, geo, &mut ctx.pad);
+ // g[i] = (f1 - f0) / step;
+ // ctx.x_tmp[i] = xx[i];
+
+ // let h = ABS_GRADIENT_STEP;
+ // ctx.x_tmp[i] = xx[i] + h;
+ // let fp = penalty(&ctx.x_tmp, s, args, geo, &mut ctx.pad);
+ // ctx.x_tmp[i] = xx[i] - h;
+ // let fm = penalty(&ctx.x_tmp, s, args, geo, &mut ctx.pad);
+ // g[i] = (fp - fm) / (2.0 * h);
+ // ctx.x_tmp[i] = xx[i];
+
+ let h = ABS_GRADIENT_STEP;
+ ctx.x_tmp[i] = xx[i] + h;
+ if CENTRAL_DIFF {
+ let fp = penalty(&ctx.x_tmp, s, args, geo, &mut ctx.pad);
+ ctx.x_tmp[i] = xx[i] - h;
+ let fm = penalty(&ctx.x_tmp, s, args, geo, &mut ctx.pad);
+ g[i] = (fp - fm) / (2.0 * h);
+ } else {
+ let f1 = penalty(&ctx.x_tmp, s, args, geo, &mut ctx.pad);
+ g[i] = (f1 - f0) / h;
+ }
+ ctx.x_tmp[i] = xx[i];
+ }
+ }
+ f0
+ };
+
+ let mut opt =
+ Nlopt::new(Algorithm::Lbfgs, x.len(), obj, Target::Minimize, ctx);
+
+ opt.set_ftol_rel(args.tolerance).ok();
+ opt.set_xtol_rel(args.tolerance).ok();
+ opt.set_maxeval(MAX_ITERATIONS as _).ok();
+
+ opt.optimize(x).ok().map(|_| x.to_vec())
+}
+fn repetition(seed: usize, args: &Args, geo: &Geometry) -> (f64, Vec) {
+ let mut rng = PcgInnerState64::oneseq_seeded(seed as u64);
+ eprintln!(
+ "Starting attempt {} with seed {}",
+ seed,
+ rng.oneseq_rxs_m_xs()
+ );
+ let mut rand = || rng.oneseq_rxs_m_xs() as f64 / u64::MAX as f64;
+
+ let mut s = (args.inner_polygons as f64).sqrt() * (2.0 + rand() * 2.0);
+ let initial_s = s;
+ let mut x = vec![0.0; args.inner_polygons * 3];
+
+ if rand() < 0.5 {
+ x.iter_mut().for_each(|v| *v = (rand() - 0.5) * s);
+ } else {
+ let gc = (args.inner_polygons as f64).sqrt().ceil() as usize;
+ let span = s * 0.45;
+ for i in 0..args.inner_polygons {
+ let step = if gc > 1 {
+ 2.0 * span / (gc - 1) as f64
+ } else {
+ 0.0
+ };
+ x[i * 3] = (i % gc) as f64 * step - span;
+ x[i * 3 + 1] = (i / gc) as f64 * step - span;
+ x[i * 3 + 2] = rand() * 2.0 * PI;
+ }
+ }
+
+ let mut bump = Bump::new();
+ let (mut best_x, mut best_s) = (x.clone(), s);
+
+ loop {
+ if let Some(mut refined) = local_min(&mut x, s, args, geo, &mut bump) {
+ let mut pad =
+ ScratchPad::new(&bump, args.inner_polygons, args.inner_sides);
+ let p = penalty(&refined, s, args, geo, &mut pad);
+ if p < args.tolerance {
+ best_x.clone_from(&refined);
+ best_s = s;
+ let base = (args.inner_polygons as f64).sqrt()
+ * args.inner_sides as f64
+ / args.container_sides as f64;
+ // let m = (s - base).mul_add(
+ // -((0.01 - args.finalstep) / (initial_s - base)),
+ // 1.0 - args.finalstep,
+ // );
+ let diff_s = s - base;
+ let range_s = initial_s - base;
+ let step_ratio = (0.01 - args.finalstep) / range_s;
+ let m = (1.0 - args.finalstep) - (diff_s * step_ratio);
+ refined.chunks_exact_mut(3).for_each(|c| {
+ c[0] *= m;
+ c[1] *= m;
+ });
+ x = refined;
+ s *= m;
+ continue;
+ }
+ if let Some(bh) = run_bh(&refined, p, s, args, geo, &mut rand) {
+ best_x.clone_from(&bh);
+ best_s = s;
+ let base = (args.inner_polygons as f64).sqrt()
+ * args.inner_sides as f64
+ / args.container_sides as f64;
+ // let m = (s - base).mul_add(
+ // -((0.01 - args.finalstep) / (initial_s - base)),
+ // 1.0 - args.finalstep,
+ // );
+ let diff_s = s - base;
+ let range_s = initial_s - base;
+ let step_ratio = (0.01 - args.finalstep) / range_s;
+ let m = (1.0 - args.finalstep) - (diff_s * step_ratio);
+ for i in 0..args.inner_polygons {
+ x[i * 3] = bh[i * 3] * m;
+ x[i * 3 + 1] = bh[i * 3 + 1] * m;
+ x[i * 3 + 2] = bh[i * 3 + 2];
+ }
+ s *= m;
+ continue;
+ }
+ }
+ break;
+ }
+ (best_s, best_x)
+}
+
+fn run_bh(
+ refined: &[f64], p: f64, s: f64, args: &Args, geo: &Geometry,
+ rand: &mut dyn FnMut() -> f64,
+) -> Option> {
+ let (mut best_p, mut best_x) = (p, refined.to_vec());
+ let (mut cur_p, mut cur_x) = (p, refined.to_vec());
+ let mut trial = cur_x.clone();
+ let mut bump = Bump::new();
+
+ for _ in 0..BH_STEPS {
+ // let trial: Vec<_> = cur_x
+ // .iter()
+ // .map(|&v| v + (rand() - 0.5) * BH_STEP_SIZE)
+ // .collect();
+ for (trial, cur_x) in trial.iter_mut().zip(&cur_x) {
+ *trial = *cur_x + (rand() - 0.5) * BH_STEP_SIZE;
+ }
+ // if let Some(opt_x) = local_min(&trial, s, args, geo, &mut bump) {
+ if let Some(opt_x) = local_min(&mut trial, s, args, geo, &mut bump) {
+ let mut pad =
+ ScratchPad::new(&bump, args.inner_polygons, args.inner_sides);
+ let opt_p = penalty(&opt_x, s, args, geo, &mut pad);
+ if opt_p < best_p {
+ best_p = opt_p;
+ best_x.clone_from(&opt_x);
+ }
+ if opt_p < cur_p || ((cur_p - opt_p) / BH_TEMP).exp() > rand() {
+ cur_p = opt_p;
+ cur_x = opt_x;
+ }
+ if best_p < args.tolerance {
+ return Some(best_x);
+ }
+ }
+ }
+ None
+}
+
+fn verify_results(v: &[f64], s: f64, args: &Args, geo: &Geometry) {
+ println!("--- Verification Pass ---");
+ let bump = Bump::new();
+ let pad = ScratchPad::new(&bump, args.inner_polygons, args.inner_sides);
+ let limit = geo.apothem * s;
+ let mut overlap_count = 0;
+ let mut out_of_bounds = 0;
+ let nsi = args.inner_sides;
+
+ for i in 0..args.inner_polygons {
+ let (px, py, pr) = (v[i * 3], v[i * 3 + 1], v[i * 3 + 2]);
+ let (sin_r, cos_r) = pr.sin_cos();
+ for j in 0..nsi {
+ let tx = cos_r * geo.inner_verts_x[j]
+ - sin_r * geo.inner_verts_y[j]
+ + px;
+ let ty = sin_r * geo.inner_verts_x[j]
+ + cos_r * geo.inner_verts_y[j]
+ + py;
+ pad.poly_verts[(i * nsi + j) * 2] = tx;
+ pad.poly_verts[(i * nsi + j) * 2 + 1] = ty;
+ }
+ for j in 0..nsi {
+ let ax = geo.inner_axes_x[j];
+ let ay = geo.inner_axes_y[j];
+ pad.poly_axes[(i * nsi + j) * 2] = cos_r * ax - sin_r * ay;
+ pad.poly_axes[(i * nsi + j) * 2 + 1] = sin_r * ax + cos_r * ay;
+ }
+ }
+
+ for i in 0..args.inner_polygons {
+ let mut is_out = false;
+ for j in 0..nsi {
+ let tx = pad.poly_verts[(i * nsi + j) * 2];
+ let ty = pad.poly_verts[(i * nsi + j) * 2 + 1];
+ if geo
+ .cont_axes_x
+ .iter()
+ .zip(&geo.cont_axes_y)
+ .any(|(&ax, &ay)| tx * ax + ty * ay > limit + VERIFY_EPS)
+ {
+ is_out = true;
+ break;
+ }
+ }
+ if is_out {
+ out_of_bounds += 1;
+ }
+
+ for j in i + 1..args.inner_polygons {
+ let mut collision = true;
+ for k in 0..(nsi * 2) {
+ let (ax, ay) = if k < nsi {
+ (
+ pad.poly_axes[(i * nsi + k) * 2],
+ pad.poly_axes[(i * nsi + k) * 2 + 1],
+ )
+ } else {
+ let k2 = k - nsi;
+ (
+ pad.poly_axes[(j * nsi + k2) * 2],
+ pad.poly_axes[(j * nsi + k2) * 2 + 1],
+ )
+ };
+ let (min1, max1) =
+ project_soa(pad.poly_verts, i * nsi * 2, nsi, ax, ay);
+ let (min2, max2) =
+ project_soa(pad.poly_verts, j * nsi * 2, nsi, ax, ay);
+ if max1.min(max2) - min1.max(min2) <= VERIFY_EPS {
+ collision = false;
+ break;
+ }
+ }
+ if collision {
+ overlap_count += 1;
+ }
+ }
+ }
+ println!(
+ "Out of bounds: {}, Pairwise overlaps: {}",
+ out_of_bounds, overlap_count
+ );
+}
+
+fn render(args: &Args, geo: &Geometry, s: f64, v: &[f64]) {
+ let path = format!(
+ "{}_{}_in_{}.png",
+ args.inner_polygons, args.inner_sides, args.container_sides
+ );
+ let root = BitMapBackend::new(&path, (1024, 1024)).into_drawing_area();
+ root.fill(&WHITE).ok();
+ let mut chart = ChartBuilder::on(&root)
+ .build_cartesian_2d(-s..s, -s..s)
+ .unwrap();
+
+ let mut c_pts: Vec<_> = (0..args.container_sides)
+ .map(|i| {
+ let a = i as f64 * (2.0 * PI / args.container_sides as f64);
+ (a.cos() * s, a.sin() * s)
+ })
+ .collect();
+ c_pts.push(c_pts[0]);
+ chart
+ .draw_series(iter::once(PathElement::new(c_pts, BLACK)))
+ .ok();
+
+ for i in 0..args.inner_polygons {
+ let (px, py, pr) = (v[i * 3], v[i * 3 + 1], v[i * 3 + 2]);
+ let (sin_r, cos_r) = pr.sin_cos();
+ let mut p_pts: Vec<_> = (0..args.inner_sides)
+ .map(|j| {
+ let tx = cos_r * geo.inner_verts_x[j]
+ - sin_r * geo.inner_verts_y[j]
+ + px;
+ let ty = sin_r * geo.inner_verts_x[j]
+ + cos_r * geo.inner_verts_y[j]
+ + py;
+ (tx, ty)
+ })
+ .collect();
+ chart
+ .draw_series(iter::once(Polygon::new(
+ p_pts.clone(),
+ RGBColor(204, 204, 204).filled(),
+ )))
+ .ok();
+ p_pts.push(p_pts[0]);
+ chart
+ .draw_series(iter::once(PathElement::new(p_pts, BLACK)))
+ .ok();
+ }
+}