[Discussion] A possible strategy for optimizing search space and symmetries

[tsrasmussen]

For an n-polycube:

All n-polycubes can be generated in a space with dimensions (n, n, n).
Most n-polycubes exist in more compact spaces.

Idea: Partition space into cells of dimension x, y, z.

Constraints:
x, y, z are all greater than or equal to 1.
Cells must be large enough to accomodate the n-polycube: x*y*z >= n
Cells should not be larger than needed: x+y+z = n+2
Permutations of x,y,z yield redundant cells.

Python function for generating unique set of cells:

from itertools import combinations_with_replacement

def find_unique_cells(n):
    """
    Finds the smallest set of cells to be populated by a polycube of size n.

    Parameters:
    n (integer): the size of the polycubes.

    Returns:
    Ordered list of tuples defining the dimensions for cells.

    """
    cells = set()
    for combo in combinations_with_replacement(range(1, n + 1), 3):
        x, y, z = sorted(combo, reverse=True)  # Sort in reverse order
        if x * y * z >= n and x + y + z == n + 2:
            cells.add((x, y, z))
    return sorted(
        cells, key=lambda group: group[0], reverse=True
    )  # Sort by the largest digit

Examples:
n = 1 | Cells: [(1, 1, 1)]
n = 2 | Cells: [(2, 1, 1)]
n = 3 | Cells: [(3, 1, 1), (2, 2, 1)]
n = 4 | Cells: [(4, 1, 1), (3, 2, 1), (2, 2, 2)]
n = 5 | Cells: [(5, 1, 1), (4, 2, 1), (3, 3, 1), (3, 2, 2)]
n = 6 | Cells: [(6, 1, 1), (5, 2, 1), (4, 3, 1), (4, 2, 2), (3, 3, 2)]
n = 7 | Cells: [(7, 1, 1), (6, 2, 1), (5, 2, 2), (5, 3, 1), (4, 3, 2), (4, 4, 1), (3, 3, 3)]
n = 8 | Cells: [(8, 1, 1), (7, 2, 1), (6, 2, 2), (6, 3, 1), (5, 3, 2), (5, 4, 1), (4, 4, 2), (4, 3, 3)]
n = 9 | Cells: [(9, 1, 1), (8, 2, 1), (7, 3, 1), (7, 2, 2), (6, 3, 2), (6, 4, 1), (5, 3, 3), (5, 4, 2), (5, 5, 1), (4, 4, 3)]
n = 10 | Cells: [(10, 1, 1), (9, 2, 1), (8, 2, 2), (8, 3, 1), (7, 3, 2), (7, 4, 1), (6, 3, 3), (6, 4, 2), (6, 5, 1), (5, 5, 2), (5, 4, 3), (4, 4, 4)]

Consider n = 4. Instead of a space of total size 4*4*4 = 64 we reduce it to (4*1*1)+(3*2*1)+(2*2*2) = 18

Cells come in 4 categories: (x=y=z), (x>y=z), (x=y>z), (x>y>z)

The Isometric or Cubic Cell: (x=y=z)
The Tetragonal Cells: (x>y=z) & (x=y>z)
The Orthorhombic Cell: (x>y>z)

Each cell category will only permit some rotations as cell dimension/shape must be preserved.
Effectively you don't want to rotate your n-polycube out of it's cell.

This should reduce the number of calculations of rotations and results to look up.

Things to consider further:

• How to ensure all unique n-polycubes are generated in a cell
• Is there a fixed way to build all?
• Could this make the crop_cube(cube) function redundant?
• Put the (n-1)-polycubes into cells for n-polycubes to generate new ones?
• Which (n-1)-polycubes go into which n-polycube cells and how?

I thought of starting at (1, 1, 1) but that will lead to some issues. with for instance not generating the 7-polycube of octahedral shape.
Should you start in the middle of the x,y-plane? In the middle of a cell?

It's also possible to generate shapes by applying symmetry operations on n-polycubes in cells. The 7-polycube of octahedral shape could be constructed from something like the (2, 2, 1) cell:
[[1 0]
[1 1]]
But I'm not sure how this would work with ensuring all possible shapes get generated.

Any comments as to why any of this shouldn't work or improve runtime and clever ideas to explore further are welcome.

[tsrasmussen]

n-polycube and n-cell inheritance scheme

n-cells can inherit (n-1)-polycubes from all (n-1)-cells that fit into n-cell.

Make a dictionary with keys of cells and values of lists of possible cells to inherit polycubes from.

cell_inheritance_dict = {}

# Put cells in dict as keys
n = 4
for i in range(1, n + 1):
    cells = find_unique_cells(i)
    for cell in cells:
        if cell not in cell_inheritance_dict:
            cell_inheritance_dict[cell] = None

# Update values for keys in dict
for key in cell_inheritance_dict:
    cell = key
    previous_cells = []
    if cell[2] > 1:
        previous_cell = (cell[0], cell[1], cell[2] - 1)
        previous_cells.append(previous_cell)
    if cell[1] > cell[2]:
        previous_cell = (cell[0], cell[1] - 1, cell[2])
        previous_cells.append(previous_cell)
    if cell[0] > cell[1]:
        previous_cell = (cell[0] - 1, cell[1], cell[2])
        previous_cells.append(previous_cell)

    cell_inheritance_dict.update({cell: previous_cells})

Example - Polycubes and cells of dimension n = 4

Current cell: (1, 1, 1) Inherits polycybes from: []
Current cell: (2, 1, 1) Inherits polycybes from: [(1, 1, 1)]
Current cell: (3, 1, 1) Inherits polycybes from: [(2, 1, 1)]
Current cell: (2, 2, 1) Inherits polycybes from: [(2, 1, 1)]
Current cell: (4, 1, 1) Inherits polycybes from: [(3, 1, 1)]
Current cell: (3, 2, 1) Inherits polycybes from: [(3, 1, 1), (2, 2, 1)]
Current cell: (2, 2, 2) Inherits polycybes from: [(2, 2, 1)]

[jamesalewis]

This was my thought as well, noting that all possible polycubes will be contained in a space roughly 1/6 the size of its bounding cube. This space, being a right tetrahedron (3D version of a right triangle) can be a 1kb datum (969 bits) for a 17-polycube. This tetrahedronal array will be n×(n+1)×(n+2)/6 cells.

We can then run it through all its possible rotations and store the lowest value, treating it like a really big int. That really big int is a lossless representation of the shape which we can use to look up in the list of found shapes.

Alternatively, we can use a hash of that datum, and store that in a lookup and double check with the source if we find a match. That hashing step may also save on memory by not needing the whole array of found shapes in memory, only their hashes, when checking for matches.

Rotating and translating to find the "minimum" will eliminate any need for rotation later.

(PS: because of treating it like a big int, any rotation we do will either totally change the int value or reveal symmetry, so we don't need to worry about having to store multiple rotational versions of a found shape.)

(PPS: In reality, that maximum cellular array can be trimmed down yet more, since a rotation of a poly cube will favor one axis over the others, and the rotationally minimum shape will never reach two, let allone all three, extremes. For example, take a 10-polycube with legs stretching out along each axis equally: this is bound not to the 10-unit tretrahedron but 4 units, and moving any of the ends further out will only find a minimum rotational position along the primary axis. I have yet do a little more thinking on this to see if we can reliably represent in a tighter array than in n×(n+1)×(n+2)/6 cells)

[jamesalewis]

Doing a little more thinking on this, we may be able to trim the cellular array further by about half. The reason is that we will find the "minimum rotation" will be contained in a space somewhat like this:

This is for a 10-polycube. With some good spatial reasoning, you'll find that any arrangement which reaches outside these bounds has a rotation that fits within it.

[MyrddinE]

If the goal is to find a unique rotation for every permutation, then it seems that this would be the appropriate algorithm:

1. Create a working space that is n^3 (naieve storage), generate the 24 permutations, and shift each permutation to the minimum corner of the cube.
2. Convert each rotation to a 1d bit array.
3. Find the 'minimum' array, when treated like an integer.

Once you've done that, you'll AUTOMATICALLY have selected the rotation that fits within the bounding box described, with the long dimension corresponding to the low (1s) coordinate, and the two medium dimensions corresponding to the 10s and 100s 'place' on the flattened coordinate.

As an example let's look at a 10 polycube. The line, when shifted to the minimum position, results in 3 unique results. The smallest is where the binary 1s extend through the 1s position, so you get {0,0,0,0,0,[985 more 0s],1,1,1,1,1,1,1,1,1,1}. The medium sized result (again, when treated like a binary number) is when the digits extend along the second axis; if the digits were numbered from the lowest value to the highest 1s would appear at position 1, 11, 21, 31, 41, 51, 61, 71, 81, and 91. In binary, it would look like {0,0,0,0,0,[905 more 0s],1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,[repeat 7 more times],1,0,0,0,0,0,0,0,0,0,1}
And finally, if it extended along the 100s axis, the digits would appear at positions 1, 101, 201, 301, etc.

The point being that if you convert all 24 rotations to binary in the n^3 cube, taking the minimum of the 24 results (viewed as a binary number) will automatically get us two things:

1. A unique rotation, so we only have to store it once.
2. A rotation that automatically fits within the minimum shape described by @jamesalewis, so we can use fewer bits to store it.

In fact, once we have determined this minimum rotation we can renumber it in binary as fitting within the minimum shape. It would need to be 'unpacked' before use (since the bit array wouldn't correspond to cuboid coordinate points).

[jamesalewis]

Another thing to note is that what I described at first is the container for all rotations of every valid polycube which intersects with the three origin planes, but the "minimum rotation" will fit in a significantly smaller space. I've done a little more thinking, and I believe a right tetrahedronal space is still the correct shape, but the three dimensions will be roughly n × ½n × ⅓n (the limit as n→∞, obviously rounding up for the finite cases) because any polycube that extends outside this space can be rotated, transposed, and/or flipped to fit within it. To demonstrate, let's consider the extremities of a 2D polysquare for n=5 (square L, Long L, and Line):

I haven't yet proven the hypothesis that all polycubes will fall within this n × ½n × ⅓n tetrahedronal space, but it appears true. This will quickly turn from a computerphile into a numberphile topic as we go into the third dimension, but I am increasingly convinced of it.

Here's the 2D cases for n=1-6, along with a persistent ID scheme (yes, I have an error in the screenshot, holdover from an ID scheme change):

EDIT: What this might look like in 3D:

[tsrasmussen]

Doing a little more thinking on this, we may be able to trim the cellular array further by about half. The reason is that we will find the "minimum rotation" will be contained in a space somewhat like this:

This is for a 10-polycube. With some good spatial reasoning, you'll find that any arrangement which reaches outside these bounds has a rotation that fits within it.

How do you construct that space for a given n? What are you boundary conditions?

Looking at the cells proposed earlier, for the n = 10 case you would get the "supercell" or space that will be guaranteed to fit all possible 10-polycubes as:

It looks as if there is an in-between with this and the other structure, that will be a bit smaller than either. For the n=10 case, it looks likes it is trimming a row of the top layer, but I'd like to know why we can do this.

The other supercell/spaces up to n = 12:

I've tried writing some code that orients the polycubes in the *.npy files in a way to force them into a rotation that prioritizes the x >= y and y >= z constraints. This effectivly shifts them to the minimum corner as MyrddinE mentioned.

The steps for a given polycube are:

1. Ensure polycube.shape : x >= y and y >= z
2. Perform all rotations and select all that have a cube in the lowest corner, polycube[0,0,0,] == 1
   If no rotation satisfies this, choose the original polycube
3. Select all rotations that still ensure polycube.shape : x >= y and y >= z
4. Select all rotations that have the highest value of ones in the first row
   • This selects all those that favour a generation scheme that prioritizes x >= y and y >= z
5. Discard identical polycubes
6. Select the polycube that has 1s at the lowest indices, indices = np.argwhere(polycube == 1)
   • This selects the polycube oriented with most cubes towards the minimum corner

I haven't thought too much of optimizing this, nor the code used which is rather messy, but was a way for me to keep track.
I've included the code below. The python render_shapes2 function is just something used to visualize the superposition of all oriented polycubes from an appropriate point of view. This superposition must fit within the supercell.

def render_shapes2(shapes, path, elev=30, azim=45, title=None):
    n = len(shapes)
    dim = max(max(a.shape) for a in shapes)
    i = int(np.ceil(np.sqrt(n)))
    voxel_dim = dim * i
    voxel_array = np.zeros((voxel_dim + i, voxel_dim + i, dim), dtype=np.byte)
    pad = 1
    for idx, shape in enumerate(shapes):
        x = (idx % i) * dim + (idx % i)
        y = (idx // i) * dim + (idx // i)
        xpad = x * pad
        ypad = y * pad
        s = shape.shape
        voxel_array[x : x + s[0], y : y + s[1], 0 : s[2]] = shape

    voxel_array = voxel_array[:voxel_dim, :voxel_dim, :]
    colors = np.empty(voxel_array.shape, dtype=object)
    colors[:] = "#FFD65DC0"

    ax = plt.figure(figsize=(20, 16), dpi=600).add_subplot(projection="3d")
    ax.view_init(elev=elev, azim=azim)  # Set the elevation and azimuth angles
    ax.voxels(voxel_array, facecolors=colors, edgecolor="k", linewidth=0.1)

    ax.set_xlim([0, voxel_array.shape[0]])
    ax.set_ylim([0, voxel_array.shape[1]])
    ax.set_zlim([0, voxel_array.shape[2]])
    ax.set_box_aspect((1, 1, voxel_array.shape[2] / voxel_array.shape[0]))

    # Add the title if provided
    if title:
        ax.set_title(title, fontsize=50)

    plt.axis("off")
    plt.savefig(path + ".png", bbox_inches="tight", pad_inches=0)

    plt.close()


def orient_polycube_for_shape(polycube):
    x, y, z = polycube.shape

    # Rotate the polycube to have x >= y and y >= z
    if x < y:
        polycube = np.rot90(
            polycube, k=1, axes=(1, 0)
        )  # Rotate 90 degrees in the x-y plane
    if y < z:
        polycube = np.rot90(
            polycube, k=1, axes=(2, 1)
        )  # Rotate 90 degrees in the y-z plane

    return polycube


def are_arrays_equal(array1, array2):
    """
    Checks if two arrays are equal element-wise.

    Args:
      array1: NumPy array.
      array2: NumPy array.

    Returns:
      True if the arrays are equal element-wise, False otherwise.
    """
    return np.array_equal(array1, array2)



def orient_polycubes(polycubes):
    best_rotations = []

    for polycube in polycubes:
        oriented_polycube = orient_polycube_for_shape(polycube)
        possible_rotations = []

        for rotation in all_rotations(oriented_polycube):
            if rotation[0, 0, 0] == 1:
                possible_rotations.append(rotation)

        if not possible_rotations:
            possible_rotations.append(oriented_polycube)

        # Filter possible_rotations based on shape constraint x >= y and y >= z
        filtered_rotations = [
            rotation
            for rotation in possible_rotations
            if rotation.shape[0] >= rotation.shape[1]
            and rotation.shape[1] >= rotation.shape[2]
        ]

        # Calculate sum of ones in the first row for each filtered rotation
        sums = [np.sum(rotation[:, 0, 0]) for rotation in filtered_rotations]

        # Find the maximum sum and corresponding rotations
        max_sum = max(sums)
        max_sum_rotations = [
            filtered_rotations[i] for i, s in enumerate(sums) if s == max_sum
        ]

        # Remove duplicates from max_sum_rotations
        unique_rotations = []
        for rotation in max_sum_rotations:
            if not any(
                are_arrays_equal(rotation, unique_rot)
                for unique_rot in unique_rotations
            ):
                unique_rotations.append(rotation)

        # Select the rotation with the most 1s at the lowest indices
        best_rotation = None
        lowest_index_sum = None
        for unique_rotation in unique_rotations:
            indices_sum = np.sum(np.argwhere(unique_rotation == 1))
            if lowest_index_sum is None or indices_sum < lowest_index_sum:
                best_rotation = unique_rotation
                lowest_index_sum = indices_sum

        best_rotations.append(best_rotation)

    return best_rotations


def generate_superposition(polycubes):
    # Find the dimensions of the bounding box that encompass all polycubes
    max_x = max(polycube.shape[0] for polycube in polycubes)
    max_y = max(polycube.shape[1] for polycube in polycubes)
    max_z = max(polycube.shape[2] for polycube in polycubes)

    # Create an empty 3D numpy array for the superposition
    superposition = np.zeros((max_x, max_y, max_z), dtype=np.byte)

    # Place each polycube in the superposition
    for polycube in polycubes:
        indices = np.argwhere(polycube == 1)
        for index in indices:
            x, y, z = index
            superposition[x, y, z] = 1

    return superposition


arrays = np.load("cubes_6.npy", allow_pickle=True)
oriented_arrays = orient_polycubes(arrays)
super_pos = [generate_superposition(oriented_arrays)]
render_shapes2(super_pos, "super_pos_6", title="Superposition of n = 6 polycubes")

It's a bit messy, but when creating a super position of the resulting polycubes I get these shapes for n = 4, 5, 6 (pardon the size of the pictures):

Aside from n = 4, the superposition of the polycubes seem to fill the entire supercell.
I'm not sure why, but I could not use the *.npy files for n larger than 6, and as of yet I have no idea why. I get the following error:

line 196, in orient_polycubes
    max_sum = max(sums)
              ^^^^^^^^^
ValueError: max() arg is an empty sequence

I'm unsure of the feasibility of a scheme implementing something like the algorithm for orienting a polycube.
You wouldn't have to check each rotation against all other polycubes, as you always choose a specific orientation, but you would have to generate all rotations of a polycube and filter those in a way that ensures uniqueness.

[MyrddinE]

@tsrasmussen, note that finding the minimum dimension (x>y>z) is not the same as finding the actual minimum orientation. It will also not guarantee a unique orientation for a given shape, making duplicate detection harder.

For example, both these (identical) n=7 polycubes are the same x, y, and z dimensions but they are not both the minimum rotation when interpreting the bit representation as an integer. Looking at your code, I did not see anything that would ensure one orientation was preferred over the other (I may have missed it) when they had the same dimensions. Without that, duplicate detection becomes more difficult.