Polyomino packer using ILP written by gemini.

packer.py

import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds


class PolyominoPacker:
    def __init__(self, grid_rows, grid_cols):
        self.R = grid_rows
        self.C = grid_cols
        self.pieces = []  # List of shapes (each shape is list of coords)

    def add_piece(self, shape_coords):
        """
        shape_coords: List of (r, c) tuples defining the shape.
        e.g., [(0,0), (1,0), (0,1)] for an L-shape
        """
        # Normalize shape to ensure top-left is (0,0)
        min_r = min(r for r, c in shape_coords)
        min_c = min(c for r, c in shape_coords)
        norm_shape = sorted([(r - min_r, c - min_c) for r, c in shape_coords])
        self.pieces.append(norm_shape)

    def _get_rotations(self, shape):
        """Generate 4 cardinal rotations of a shape."""
        rotations = []
        current = shape
        for _ in range(4):
            # Normalize current rotation
            min_r = min(r for r, c in current)
            min_c = min(c for r, c in current)
            norm = sorted([(r - min_r, c - min_c) for r, c in current])

            if norm not in rotations:
                rotations.append(tuple(norm))

            # Rotate 90 degrees clockwise: (r, c) -> (c, -r)
            current = [(c, -r) for r, c in current]

        return set(rotations)

    def solve(self):
        # 1. Generate all valid placements (columns of our matrix)
        # placement_map stores metadata for each column index: (piece_idx, cells_occupied)
        placements = []

        print(
            f"Generating placements for {len(self.pieces)} pieces on {self.R}x{self.C} grid..."
        )

        for p_idx, shape in enumerate(self.pieces):
            unique_rots = self._get_rotations(shape)

            for rot_shape in unique_rots:
                max_r = max(r for r, c in rot_shape)
                max_c = max(c for r, c in rot_shape)

                # Try every translation
                for r in range(self.R - max_r):
                    for c in range(self.C - max_c):
                        # Calculate global occupied cells
                        occupied = []
                        for dr, dc in rot_shape:
                            occupied.append((r + dr, c + dc))

                        placements.append(
                            {"piece_id": p_idx, "cells": occupied}
                        )

        num_vars = len(placements)
        print(f"Total potential placements (variables): {num_vars}")

        # 2. Build Constraints

        # --- Constraint A: Each piece placed exactly once ---
        # Rows = Number of pieces
        A_pieces = np.zeros((len(self.pieces), num_vars))
        for j, p in enumerate(placements):
            A_pieces[p["piece_id"], j] = 1

        lb_pieces = np.ones(len(self.pieces))
        ub_pieces = np.ones(len(self.pieces))

        # --- Constraint B: No Overlaps ---
        # Rows = Number of grid cells
        A_grid = np.zeros((self.R * self.C, num_vars))
        for j, p in enumerate(placements):
            for r, c in p["cells"]:
                grid_idx = r * self.C + c
                A_grid[grid_idx, j] = 1

        lb_grid = np.full(self.R * self.C, -np.inf)
        ub_grid = np.ones(self.R * self.C)  # Max 1 per cell

        # Combine
        A = np.vstack([A_pieces, A_grid])
        lb = np.concatenate([lb_pieces, lb_grid])
        ub = np.concatenate([ub_pieces, ub_grid])

        # 3. Solve
        c = np.zeros(num_vars)  # No cost, just feasibility
        res = milp(
            c=c,
            constraints=LinearConstraint(A, lb, ub),
            integrality=np.ones(num_vars),
            bounds=Bounds(0, 1),
        )

        if res.success:
            print("Solution found!")
            self._visualize(res.x, placements)
        else:
            print("No valid configuration found.")

    def _visualize(self, x, placements):
        grid = np.full((self.R, self.C), ".", dtype=object)
        selected_indices = np.where(x > 0.5)[0]

        # Symbols for pieces: 0=A, 1=B, 2=C...
        symbols = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"

        for idx in selected_indices:
            p = placements[idx]
            symbol = symbols[p["piece_id"] % len(symbols)]
            for r, c in p["cells"]:
                grid[r, c] = symbol

        # Print Grid
        print("\nGrid State:")
        for row in grid:
            print(" ".join(row))


if __name__ == "__main__":

    solver = PolyominoPacker(8, 8)

    shapes = {
        "t": [(0, 1), (1, 0), (1, 1), (1, 2)],
        "l": [(0, 0), (1, 0), (2, 0), (2, 1)],
        "sq": [(0, 0), (0, 1), (1, 0), (1, 1)],
        "line": [(0, 0), (0, 1), (0, 2), (0, 3)],
        "s": [(0, 0), (1, 0), (1, 1), (2, 1)],
    }

    # How many of each shape
    shape_counts = {
        "t": 4,
        "l": 3,
        "sq": 2,
        "s": 4,
        "line": 2
    }

    # 3. Add arbitrary number of arbitrary shapes
    for key in shape_counts:
        for _ in range(shape_counts[key]):
            solver.add_piece(shapes[key])

    # 4. Run
    solver.solve()