mikepound/opencubes

A community improved version of the polycubes project!

Polycubes

This code is associated with the Computerphile video on generating polycubes. The original repository may be found here. That version is unchanged from my original video, so that those watching for the first time can find and use the original code, and make improvements to it themselves. This repository is for those looking to contribute to a faster and better optimised version, driven by improvements from Computerphile viewers!

Introduction

A polycube is a set of cubes in any configuration in which all cubes are orthogonally connected - share a face. This code calculates all the variations of 3D polycubes for any size (time permitting!).

Contents

This repository contains 3 solutions written in 3 languages Python, C++, and Rust. each sub folder contains a README with instructions on how to run them.

Improving the code

This repo already has some improvements included, and will happily accept more via pull request. Some things you might think about:

• The main limiting factor at this time seems to be memory usage, at n=14+ you need hundereds of GB's just to store the cubes, so keeping them all in main memory gets dificult.
• Distributing the computation across many systems would allow us to scale horizontally rather than vertically, but it opens questions of how to do so without each system having a full copy of all the cubes, and how to manage the large quantities of data.
• Calculating 24 rotations of a cube is slow, the only way to avoid this would be to come up with some rotationally invariant way of comparing cubes. I've not thought of one yet!

Contributing!

This version welcomes contributors!

References

• Wikipedia article
• This repository was a source of inspiration, and a great description of some possible ways to solve this.
• There may be better ways to count these, but I've not explored in much detail.
• Kevin Gong's webpage on enumerating all shapes up to n=16.

opencubes's Issues

OpenCubes file formats

There are few file formats around. It is however bit unclear what these formats actually are for.
This issue is for discussion of the file formats used project wide + their documentation.

• What is the current export file format?
  I.e. the format all implementations are expected to be able write that contains the polycubes of N=?
• What are the currently used (other) file formats of the implementations?

@bertie2 @datdenkikniet @nsch0e @NailLegProcessorDivide

Rust and python versions seem to implement an format that follows #17 (comment) so I think this is the export file format.

Rust implementation's "canonical" isn't "pcube canonical"

See title.

Keeping this open as a tracking issue, fixed in #21

Optimizations from guy that did n=16

Here is the paper from the guy calculating to n=16 in 2005 in 11 days: http://kevingong.com/Polyominoes/ParallelPoly.html

Something very interesting he wrote about how to keep the hashtables small and more dividable

4.1.2 Hashing

One of most fundamental tasks of the polyomino program is to support uniqueness. We have to ensure that a polyomino has not already been counted. To do this, we store all of the unique polyominoes in a hash tree/table structure, and compare the candidate polyomino with only a small subset of the unique polyominoes. To access the hash structure, we first branch on some criteria(on) to find the correct hash table, then index into the hash table.
If we compare the candiate polyomino only with polyominoes of the same height and width, this narrows down the search considerably. We can also branch on the y-coordinate of the last square. This seems a bit strange, but it works well. (See Figure 4.1.2-1) To narrow it down even further, we use a hash function. This takes the XOR of every n bits, where n = height-1. (See Figure 4.1.2-2) This is not necessarily the best hash function, but it is simple to implement, and gives fair results, as we will see.

Shape encoding

I had an idea for a rotation, reflection invariant method of encoding. I don't have a lot of time to try and implement it right now but I thought I would put it here. It may be expensive in itself but If it reduces the algorithmic complexity it would be worth it.

A single binary 3d matrix can have every possible transformation applied to it so that we have a set of every orientation of the shape as matrices. That set is unique to that shape. If each entry is then flattened, then sorted, and then concatenated you will have a big binary lump that's unique. This can then have a standard hashing function applied to it. It could probably be compressed considerably as well.

Could we replace speed with a distributed system?

Instead of trying to compact the memory layout of the cubes is it possible to use a shared database to store the hashes of each shape? This almost feels like a blockchain style problem where everyone tries to guess the next hash, in this case it's just the next shape.

We could even make this a big distributed Numberphile viewer calculation!

Idea: New way to encode shapes?

Encoding shape by "flattening" all three dimensions

The goal is to encode shapes by having a sort of "map" that tells you how many blocks are at a certain coordinate in a certain dimension.

An example will hopefully make it clearer:

Imagine the well-known T-Piece (4 Cubes) in a specific orientation, and the "Views" from each dimension:

Base

Top view

Side view

Front view

Note the numbers i have added to the images, these denote how many cubes are behind the one face that is visible when viewing from the indicated side.

These are the essence of my proposal, so the T-Piece could be represented by the following Matrices

[1 2 1]

[0 1 0]
[1 1 1]

[1]
[3]

As far as I can tell, these three arrays will describe any orientation of the T-Piece when viewed from the correct side and you won't have to check all possible 24 rotations to compare two polycubes.

Now for my questions:

Is this feasible to implement? Will this even be more efficient than what is done right now?

Please ask any questions you might have, I'm willing to concede if turns out to be a bad idea

Representations for n=1 to n=4

n=1

Just a 1x1x1 cube
[1]
[1]
[1]

n=2

1x1x2
[1 1]
[1 1]
[2]

n=3

Corner	Long
[1 1] [0 1]	[1 1 1]
[2 1]	[1 1 1]
[1 2]	[3]

n=4

Square	L	Long	Z	T	Corner/Pyramid	"Twisted Z" (the one that can be reflected)
[1 1] [1 1]	[1 1 1] [0 0 1]	[1 1 1 1]	[1 1 0] [0 1 1]	[0 1 0] [1 1 1]	[1 2] [0 1]	[2 1] [0 1]
[2 2]	[3 1]	[1 1 1 1]	[2 2]	[3 1]	[1 0] [2 1]	[1 0] [2 1]
[2 2]	[1 1 2]	[4]	[1 2 1]	[1 2 1]	[0 1] [1 2]	[1 0] [1 2]

This should prove that each unique shape has three arrays which make a unique set to describe that shape, at least for n <= 4, hopefully for other n as well.

My rotationally invariant hashcode solution. Hash Collision Rate of 0.234% at N = 13.

As the title says, I have created a hashcode which I think is pretty good. Its not only rotationally invariant but also reflectionally invariant too, so the same polycube given in any rotation/reflection always gets the same hashcode. The dream when I set out was to make a perfect hashcode with no collision, but I just couldn't find enough rotationally invariant distinguishing data to tell all polycubes apart.

The hashcode works by calculating each cubes manhattan distance (gird moves) and euclidean distance (straight line) to every other cube, and combining them all together to make a single hash. Thats still an N squared operation, but if it avoids having to do the 24 rotations its still faster for all N's less than 24. This solution was able to achieve N = 13 in 9 minutes on my machine.

https://github.com/BenjaminFausett/PolycubeEnumerator

No hash table needed

Here's an algorithm that is able to identify unique cubes without using a hash table:

1. Start with a polycube p
2. For each possible expansion q of p:
   1. Compute all possible ways of removing one cube from q and have it still be connected, call this set S
   2. If p has the minimal canonical representation out of all polycubes in S, save q

We already have steps 1 and 2. I'm not sure how expensive steps 2 i. will be, but this should be parallelizable/distributable as well as avoiding any memory issues. Could be worth playing around with - I won't have time this weekend but I may try this out the weekend after if somebody else doesn't want to try this first.

Representation and orientation by prime products

Here's a proposal that I hope to eventually start implementing, but feel free to run with the idea if you wish.

This concerns the core topics of how to efficiently look up whether a polycube candidate is a rotation (including the null rotation) of any in the set of canonical polycubes, and how to possibly iterate over candidates in a way that lets us skip iterations/branches early if it can be concluded that they would not lead to valid or canonical polycubes.

One way to define the canonical orientation of a polycube (once it has been cropped to a tight bounding block) is to pick the orientation that gives the minimum hash value. If we can construct the hashing function to have some nice properties, it might be possible to know that candidates of size N are iterated in a canonical order, predict whether an cube-extension would give a non-canonical orientation and then start taking shortcuts in the iteration...

I propose the idea to assign a prime number to each 3D array element and then and define the hash number as the product of the primes whose locations are occupied with cubes. I.e. if occupied is a 3D-array with 1 where cubes are located and 0 elsewhere and primes is another 3D-array, the prime product would be P = np.product(primes ** occupied).

Consequences for iterating

If we for the moment ignore the issue of how many and large numbers would be needed, some nice properties for the iteration appear from the definition of canonical orientation as having the smallest P. I defined the order by defining primes = primes_1D[:H*W*L].reshape(H, W, L) for the case of a search space of shape (H, W, L) and with primes_1D being the sequence of enough primes, e.g. [2, 3, 5, 7, 11...]. Since this means that primes[0, 1, :] > primes[0, 0, -1] and so on, the canonical orientation keeps the occupied bounding box as close to the origin as possible (no padding before) and prefers to extend along along the last axis, i.e. keeping h <= w <= l for the bounding block of height h, width w and length l. Some pseudocode made me convinced that there are ways to skip many polycube extension candidates just by comparing the coordinates and whether the h <= w <= l invariant would be invalidated.

When extending towards positive coordinates or filling internal voids, I believe the canonical new P_candidate can be computed as simply as P_smaller_polycube * primes[where the new cube is put] without having to do any array padding or rotation! For extending towards negative coordinates, a special case per axis can be used, and here it seems the various rotations need to be computed to guarantee finding the canonical rotation of the candidate.

I don't have a very strong proof that my iteration scheme will find all polycubes of size N when starting from only canonically oriented cubes of size N-1, this remains to be shown empirially by comparison with the results of others...

Consequences for hashing and equality check

If we never reuse primes within the 3D array, the factorization of a product P is unique so the hashing is actually a bijective function between canonical polycube and its corresponding integer. Thus there's no need for storing something like a RLE-encoded polycube in memory to check for equality in the case of hash-collisions, we just need to count how many unique hash numbers we find, while ensuring we only do it for polycubes in canonical orientation. We can of course still write (append) newly found polycubes to a cache file in some format.

Size of numbers

Whether this would be faster than existing approaches is of course not obvious. If the prime product can fit in an unsigned 64-bit integer I have some hope. But so far I have not been able to prove that less than 173 bits would suffice for N=17 (but N=9 should definitively fit). I know my estimate is an unnecessarily high bound and Python does have arbitrary-size integers that could be used in a first implementation for exploring the range of canonical P-values empirically (letting the hash map use some internal non-unique machine-sized integers). I did consider optimizations like leaving zeros in some far-off elements of primes_3D known to not be needed for canonically oriented polycubes, and even special-case treatment of very straight polycubes to let the 3D-array side length be at most N-2. (Hash codes that involve squares of prime numbers are available for special meanings, since occupied.max() being 1 means that squares will never appear in P!) But to ensure my bit-bound is not an underestimate I currently allow prime-products that are clearly unreachable for canonically oriented cubes.

The day after first posting this, I edit to include a minor optimization towards smaller products: Instead of just using primes as factors, we can allow prime**2, prime**4, prime**8 and so on (powers of powers of 2) since there's still no way for two polycubes to get the same product. instead of primes = [2, 3, 5, ...] we now use numbers = [2, 3, 4, 5, 7, 9, 11, 13, 16, 17, ...] and for instance only the straight polycube of length 3 will have the product P = 24 = 2**1 * 3**1 * 2**2. Put differently: if we previously needed to use the V = len(primes) smallest primes, we can instead use the V smallest numbers in the set of primes**(2**k) for k < np.log2(primes.max()). For N=17, this just improves my upper bound from 173.4 to 171.2 bits, but even a two bit reduction might for some N make the difference between fitting in 64 bits or not.

the solution is:

using the "hilbert curve" on the code, then change the curve to the line, then with "catalan equation" check validation patterns.
this problem fixed years ago with mathematical proof.
cya

for help to understand can to watch this 3 min on youtube: https://www.youtube.com/watch?v=kaInaIUABzY

if after 5 steps can not growth on catalan number, use this idea:
fact check: Levine Sequence "sum of each row show the Catalan number"
Levine Sequence : https://oeis.org/A011784

Chat-GPT's Solution To A Unique Hashset Key For Each Shape (Regardless of Reflection & Rotation)

Solution

import numpy as np

def normalize_shape(shape):
    # Convert the shape to a binary numpy array
    binary_shape = np.array(shape, dtype=bool)

    # Generate canonical representations for rotations and reflections
    canonical_representations = []
    normalized_shape = binary_shape

    for _ in range(4):
        normalized_shape = np.rot90(normalized_shape)
        canonical_representations.append(normalized_shape)

    normalized_shape = np.fliplr(binary_shape)
    for _ in range(4):
        normalized_shape = np.rot90(normalized_shape)
        canonical_representations.append(normalized_shape)

    # Find the lexicographically smallest representation
    canonical_representations.sort(key=lambda x: tuple(x.flatten()))
    return canonical_representations[0]

def calculate_distance_signature(shape):
    # Convert the shape to a binary numpy array
    binary_shape = np.array(shape, dtype=bool)

    # Calculate the centroid (center of mass)
    rows, cols = np.nonzero(binary_shape)
    centroid = np.mean(rows), np.mean(cols)

    # Calculate the distances from each boundary point to the centroid
    distances = set()
    for r, c in zip(rows, cols):
        distance = np.sqrt((r - centroid[0])**2 + (c - centroid[1])**2)
        distances.add(distance)

    # Sort the distances for a consistent hash key
    sorted_distances = tuple(sorted(distances))

    return sorted_distances

def is_shape_stored(shape, hashmap):
    # Get the normalized representation of the shape
    normalized_shape = normalize_shape(shape)

    # Get the distance signature for the normalized shape
    distance_signature = calculate_distance_signature(normalized_shape)

    # Check if the distance signature is already stored
    if distance_signature in hashmap:
        return True

    # If the distance signature is not in the hashmap, store it
    hashmap.add(distance_signature)
    return False

stored_shapes = set()

# Test cases
m = [[1, 1, 0],
     [0, 1, 0],
     [0, 1, 0]]

n = [[0, 0, 0],
     [0, 0, 1],
     [1, 1, 1]]

print(is_shape_stored(m, stored_shapes))  # Output: True
print(is_shape_stored(n, stored_shapes))  # Output: True

Explanation

Problem Statement:
The problem is to determine if two 2D shapes are the same, regardless of rotation and reflection.

Solution Approach:
To achieve this, we need to find a way to convert each shape into a unique hash key such that equivalent shapes will produce the same hash key. This unique hash key allows us to perform quick lookups to check if a shape has already been encountered.

Step 1: Normalize Shapes
To handle rotations and reflections, we need to normalize the shapes into their canonical representations. We achieve this by generating all possible rotations and reflections of the shape and choosing the lexicographically smallest representation among them.

Step 2: Calculate Distance Signature
The next step is to compute the distance signature of the normalized shape. The distance signature represents the shape as a set of distances from each boundary point to the centroid (center of mass) of the shape. By using distances, we can capture the overall shape geometry and its invariance under rotation and reflection.

Step 3: Calculate Centroid and Distances
The centroid of a shape is computed as the average of the row and column indices of all the True elements (boundary points) in the shape. This gives us the center point around which we calculate the distances.

For each boundary point (row, col) in the shape, we calculate the Euclidean distance from that point to the centroid. These distances are collected in a set to eliminate duplicates.

Step 4: Sort Distances
To ensure a consistent hash key, we sort the distances obtained in Step 3. This is important because the distances can be in different orders for different shapes, but equivalent shapes should have the same distances when sorted.

Step 5: Using Hashsets for O(1) Lookup
We use a Python set (hashset) data structure to store the calculated distance signatures. Checking for the existence of a distance signature in the set has an average-case time complexity of O(1). This allows us to quickly determine if a shape has been encountered before.

Step 6: Check for Shape Existence
In the is_shape_stored function, we take a shape as input. First, we normalize the shape to get its canonical representation. Then, we calculate the distance signature of the normalized shape.

Step 7: Check for Existence in Hashset
We check if the calculated distance signature is already present in the hashset. If it is, this means the shape has been encountered before and is equivalent to a previously seen shape. Hence, we return True. Otherwise, we add the distance signature to the hashset and return False.

Conclusion:
The solution uses the normalized shape's distance signature as a unique hash key to perform shape existence checks in a hashset. This approach allows us to quickly determine if two shapes are the same regardless of their rotation or reflection. While the time complexity of calculating the distance signature is O(N) (N being the number of non-zero elements in the shape), the overall average-case time complexity of the solution remains O(1) due to the use of hashset for quick lookups.

Suggestion for cannonical representation

Hello All,

I am a amateur mathmatician who got interested by the puzzle.

What we are looking for is a way to generate all extensions of the previous unique shapes required to have generated all unique shapes of the next size.
OR
To generate all combinations of cubes required to have generated all unique shapes of a given size.

The hardest part of the problem is having to do a "contains" operation on the set of all currently found unique shapes, as it requires an expensive comparison N times, where N is the amount of shapes already found.

It would be nice if we could have a relatively cheap operation that turned a shape into it's canonical form and put it into an ordered table/hashtable, which would be an existence table/frequency table initialized at 0. saving us from ~N^N expensive operations.

We need a measure independent of rotation and translation, as all we care about is shape.

So, what if we could find a cube as a starting point in a way based on shape, rather than specifics, and then use that cube as a root node and then generate a tree structure as if it would be generated from the root cube?

If this tree structure can always be generated and always returns the same canonical form, it is a useful canonical function.

Selecting this cube might be done by finding the closest cube to the "center of mass"+ a function to generate "Up" and "Left" to be able to create an unambiguous tree structure.

Up and Left could be generated from the center of mass stat if it is biased from the true center of a cube.
If the center of mass is on a center plane of a cube, you need another exemption or full check.
If the center of mass is dead center, you need a full check.

This solution can be executed perfectly parallel.

My question to the more experienced:
How expensive is a center of mass calculation? If you work of previous data, you might you be able to fudge mass calculations.

What have I misunderstood about the complexity of the problem?
Because it feels like my suggestion is different from others, so i must have misunderstood something. Maybe center of mass is really expensive.

I look forward to your replies.

Enumerating 3d and 4d rotations at the same time

It's possible to count the 3d rotation polycubes when generating the 4d rotation polycubes using the following relationship:

If a 4d rotation polycube has a reflectional symmetry, then it counts as one 3d rotation polycube, else it counts as two.

I have a repo that implements this idea, #11, #49, and a streaming articulation point algorithm. The performance seems competitive, but I would appreciate if someone else could benchmark it as well. For n = 12 my laptop takes 5.6s to run mine, 7s to run #49 with ./polycube_generator 12 -a, 11.2s to run #27, and 39.6s to run the current rust version with cargo run -r -- enumerate -m hashless 12 and a n = 11 cache file.

Cuboids

Hello there,

This past week I have been trying to tackle a certain problem : the computation of polycubes
that fits inside a 3x3x3 cuboid. My work brought me to improve the Python code to make the computation faster
but it also gave me an idea for the computation of polycubes of size n.

I consider as valid polycubes inside a cuboid of size (k1,k2,k3),
those which span across the full extent of the cuboid in the 3 directions.
I also put the constraint of k1>=k2>=k3 which does not change the final result because the problem is rotation invariant.
The direct consequence is that the number of cubes (let me call that the rank)
n of a valid polycube is bounded like so:
k1+k2+k3-2 <= n <= k1k2k3. For the (3,3,3) cuboid, we have 7 <= n <= 27.

The generation of polycubes of each increasing rank inside a cuboid of a given shape is done starting from the rank 1 polycube
at all possible positions inside the cuboid. The expansion function is modified since the output shape is the same as the input shape.
Restricting the polycubes that fits inside the cuboid makes for a few optimizations in the computation:
First, the shape is fixed so it does not need to be added to the packing of the cube.
With Python, I have used the tobytes() function for an even faster packing but that does take 8 times as much space as needed.
Second, the number of rotations in a cuboid is limited to 8 when k1=/=k2=/=k3 which is the most occuring general situation.
In the case of the (3,3,3) cuboid, this is still 24 though.
Third, the shape being fixed, a rotation is an index permutation operation and is made much faster (at least with numpy).

As a result, the number of polycubes of each rank tends to increase exponentially from 1 till n=k1k2k3/2 and then decreases exponentially to 1.
The screening of invalid polycubes removes around 10% of the total after the expansion is over.
By the way, the result for cuboids with k3>=2 and k1k2k3<=28 are available on my online viewer.
https://asliceofcuriosity.fr/assets/polycube-viewer/index.html
I'll make a PR soon with these improvements.

As for my contribution to the computation of polycubes of rank n, here is my idea.
When you have all the polycubes in a cuboid of a certain shape (k1,k2,k3), you can compute part of the polycubes
of the cuboids of shape (k1+1,k2,k3), (k1,k2+1,k3), (k1,k2,k3+1).
Starting with the (1,1,1) cuboid, this makes for a graph of dependency where each cuboid
can be computed by the expansion in different directions of three previous cuboids
and can be expanded into three cuboids.
So if one wants to generate and enumerate all of the polycubes of rank n,
it needs to gather the polycubes of rank n from cuboids that satisfy the bounding equation
and all of the cuboids before them in the graph up to a certain rank and so on.
This dependency of cuboids and ranks is "easy" to compute and would split up the computation
into several smaller chunks which would reduce file sizes, RAM usage and maybe computation time too.

I think that grouping together polycubes with the same shape would allow for many optimizations in terms of storage and computation.
I have seen the idea mentioned in different ways across the posts from different people on different topics.
I am in the process of implementing this idea in Python.
The graph starts like this, with of all cuboids of the same depth in the graph on the same line and with the shape, the minimum and maximum rank for each cuboid:
(1, 1, 1) 1 1
(2, 1, 1) 2 2
(3, 1, 1) 3 3 / (2, 2, 1) 3 4
(3, 2, 1) 4 6 / (4, 1, 1) 4 4 / (2, 2, 2) 4 8
(5, 1, 1) 5 5 / (3, 3, 1) 5 9 / (3, 2, 2) 5 12 / (4, 2, 1) 5 8
(6, 1, 1) 6 6 / (5, 2, 1) 6 10 / (4, 3, 1) 6 12 / (3, 3, 2) 6 18 / (4, 2, 2) 6 16
(4, 3, 2) 7 24 / (5, 2, 2) 7 20 / (7, 1, 1) 7 7 / (4, 4, 1) 7 16 / (3, 3, 3) 7 27 / (6, 2, 1) 7 12 / (5, 3, 1) 7 15

Standard cube format.

Given we now have a second implementation language tentatively waiting as a pull request. it would be optimal to agree on a standard format for cubes at some point so that multiple implementations can share data.
could any ideas for this format go here please.

C solver with new dictionary key generation method

I've created a solver that uses a dictionary with keys based on rotation independent features of the cube shapes.

While the keys don't uniquely identify the shapes, they limit the scope of the search to the point where it can find the 10 cubes shapes in 47 seconds.

if you are interested it is here: https://github.com/tomestopen/CubeShapes/tree/1.0

There is definitely room for improvement, especially in the multithreaded version of the search function.

Making a distributed system, n=18 and beyond

I know @bertie2 and @datdenkikniet have also both also mentioned the possibility of a distributed system which the current DFS solutions should theoretically lend them self to extremely well.

This seems like something that comes with a lot of additional problems and would need some coordination between people to develop a concrete api as well as if we want to have at least basic guarantees for data integrity or security.

What I'm mainly asking is how would such a solution be arranged and does anyone have any experience with similar systems?

On the use of sparse Matrices

i believe using sparse matrices rather than using 3D matrices should reduce the amount of stored data by a fair amount.

No hash table with less performance penalty

Starting with the "hashtable-less" algorithm described by presseyt in #11, it seems it could be faster if we change the approach in checking if we have the correct polycube for output:

1. Start with seed polycube p. Extend by candidate cube a to yield polycube q (q = p + a).
2. Repeat step 1 for all possible additions of a.
3. Filter the results of step 2 to only have unique values, saving the highest found value of the index of a in each q.
4. For each q, check if there is any possible r at a higher index of q that can be removed while remaining a polycube. Otherwise output q.

The advantage of this method is we don't need to recanonize the resulting polycubes after removing the values of r. This doesn't come completely for free, as we have to check for duplicates within the set of polycubes generated by the seed polycube p and do some connection checks. This is because symmetric seed polycubes will create at least 2 of the same generated polycube, though we check only within this set because the index of the added point is always the same in the output for the same seed polycube.

I made my own C implementation based on this idea in this repository. Profiling on low numbers of N seems to indicate the connection checks for step 4 take about 15% of the total processing time. Testing shows it seems faster (N=15 in 1 hour 38 minutes on an old FX-8350 processor) than the current rust hashless or point-list versions on my machine, though I'm not sure how it compares to the solution described in #27. Admittedly, I'm rather limited in testing larger values of N because of how long I'll have to tie up my machine with its limited processing capability.

Hierarchical pcube representation

Following is just me thinking about the polycube representation and how an pcube is computed.
(sorry if the grammar is bad/gibberish as I'm not too good writing complex sentences in english... 🥲 )

An tree (hierarchical representation) is an compact way to represent the polycubes:

When an base pcube N is expanded to N+1 etc. we count the index of the permutation:

1. Generate unique list of candidate expansions for the base pcube.
2. Start with single candidate in the permutation list.
   This is the "normal way" N+1 expansion of pcube.
3. Process all unique permutations of the expansion list.
   For N+1 expansion the permutation index is just the candidate's array index.

Whats new is the N+2, N+3, etc. expansion of an base pcube:

• N+2: Generate all unique permutations with two candidates.
  Each unique shuffle increments the permutation index.
• N+3 expansion is all unique permutations with three candidates.
• ...
• The final base pcube expansion is simply using all the candidates.

Above process gives unique "index" value to all possible base pcube expansions:
Since we are computing just N+1 the expansions beyond this are not needed until later.
(Note: if we remember the last permutation state for base pcube N+2,
we can continue the enumeration from where we left it for future N+3 etc.)

Now that base pcube expansions have an well defined "index" encode pcube(s) with following:

• Some kind of reference to the base pcube
• the used permutation index value.
  Note that this value will grow with the count of candidate combinations.
• the rotation's index. (an number of [0, 23] )

This scheme does mean that to get the actual coordinate data we have to (re-)compute it using
the expansion process again:

1. Start with already rotated base pcube.
2. Compute expansion candidate list for this pcube.
3. Expand using the permutation index.
   This an hard problem: If we have to start from zero this would take an long time.
   Better options would be to selectively save the permutation state+index and continue from nearest check point.
4. Rotate the expanded pcube using the rotation index.
5. Repeat until we have pcube of required N

The hierarchical representation can be collapsed this way into the canonical form. (sorted coordinate list)

The reference to the base pcube de-duplicates it for all pcubes derived from it and this forms an tree structure.
The trick for saving memory is that the base pcube(s) should be loaded from the disk cache on demand
so N-1 part of the pcube(s) is kept out of memory.
Draw back is that the disk cache system needs to support finding any particular pcube and we have to store the pcube tree structure on disk.

Does any above idea(s) make sense ?

Rotation is a significant factor in runtime

From profiling the function generate_polycubes, for larger values of n the time spent on rot90 is a significant factor, with 45 of the 108 second runtime being spent in it for n=9 (see table at bottom of text).

So far I've been thinking about a couple of potential solutions, but I haven't reach a point where either solution could be implemented:

Solution A - Writing a purpose-specific rot90 or all_rotations

Numpy appears to use index swapping to achieve 90 degree rotation of an array (definitely faster than a transformation matrix I would imagine), but it has significant overhead as it is written to work for n-dimensional arrays. In the case of polycubes, we know that the array we want to rotate is 3D, and that we want all 28 rotations. Writing a purpose specific version of this function could save on overhead time, but I think it would most likely need to be done in another language.

I haven't really been pursuing this as I don't know any low level or performant languages to try this in, and I think we won't save time with a pure python implementation.

Solution B - Representing the polycube as three 2-D views

This solution appeals to me far more than solution A, as it reduces either our comparisons between hashes, or our memory used on storing hashes (if we store rotations). As we conveniently have a boolean array aligned with the axes, I think we can represent any polycube with three, fast to calculate, 2-D views of the polycube:

def view_2d(polycube):
    view1 = np.sum(polycube, axis=0)
    view2 = np.sum(polycube, axis=1)
    view3 = np.sum(polycube, axis=2)
    
    return view1, view2, view3

For the simple L shaped (cropped) polycube with a length of 4 and 1 cube off the a side, the views would be represented this way:

[4, 1], [1, 1, 1, 2], [1, 0, 1, 0, 1, 0, 1, 1]

We can avoid worrying about the order in which these appear by using collections.Counter or generally by treating it as an unsorted list. We could then compare the view to a rotated version of the polycube by ensuring that all entries match only once, hence finding that both polycubes are identical...

Unfortunately, though this works for comparing the two 1-D views of a polysquare with some additional checks, I haven't found a way to ensure these views match for polycubes after any of the 28 arbitrary rotations:

polycube1_views         -> [4, 1], [1, 1, 1, 2], [1, 0, 1, 0, 1, 0, 1, 1]
polycube1_rotated_views -> [1, 4], [1, 0, 1, 0, 1, 1, 1], [1, 1, 1, 2]

We can see from this that some rotations will give us outputs that are very hard to compare to the original polycube for equality.

Perhaps a solution to this is to rebuild the arbitrary polycube in a specific manner to always give the same output? I don't have a solution for this yet.

Profiling

cumtime is the value we're worried about for rot90 in this case, as looking at the source code for this function, it calls a number of other functions to achieve rotation by 90 degrees, with flip being the main point of concern.

With optimisations from the PRs on the original repo (bertie2, georgedorn, and mine), the cProfile.run output looks like this for n=9:

         67515198 function calls (64901683 primitive calls) in 108.627 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   669114    1.084    0.000    8.019    0.000 <__array_function__ internals>:177(any)     
     8151    0.013    0.000    0.055    0.000 <__array_function__ internals>:177(around)  
        1    0.000    0.000    0.000    0.000 <__array_function__ internals>:177(copyto)  
  1693830    2.401    0.000   21.180    0.000 <__array_function__ internals>:177(flip)    
  1578102    2.321    0.000    4.637    0.000 <__array_function__ internals>:177(packbits)
     8151    0.018    0.000    0.794    0.000 <__array_function__ internals>:177(pad)     
  1637369    2.486    0.000   45.442    0.000 <__array_function__ internals>:177(rot90)   
     8151    0.014    0.000    0.091    0.000 <__array_function__ internals>:177(round_)  
  1338228    1.863    0.000    6.435    0.000 <__array_function__ internals>:177(swapaxes)
   903376    1.308    0.000    6.317    0.000 <__array_function__ internals>:177(transpose)
        1    0.055    0.055  108.627  108.627 <string>:1(<module>)
    32604    0.023    0.000    0.023    0.000 arraypad.py:109(<genexpr>)
    32604    0.021    0.000    0.021    0.000 arraypad.py:120(<genexpr>)
    24453    0.125    0.000    0.155    0.000 arraypad.py:129(_set_pad_area)
    48906    0.030    0.000    0.030    0.000 arraypad.py:33(_slice_at_axis)
    16302    0.086    0.000    0.220    0.000 arraypad.py:454(_as_pairs)
     8151    0.002    0.000    0.002    0.000 arraypad.py:521(_pad_dispatcher)
     8151    0.169    0.000    0.746    0.000 arraypad.py:529(pad)
    24453    0.033    0.000    0.033    0.000 arraypad.py:58(_view_roi)
     8151    0.105    0.000    0.157    0.000 arraypad.py:86(_pad_simple)
  1578102    8.903    0.000   17.019    0.000 cubes.py:152(rle)
   223038    0.341    0.000    2.796    0.000 cubes.py:174(cube_exists_rle)
  1693830    4.488    0.000   40.548    0.000 cubes.py:24(single_axis_rotation)
   223038    7.844    0.000   22.298    0.000 cubes.py:44(crop_cube)
   231189    1.430    0.000   24.981    0.000 cubes.py:64(expand_cube)
  1411525    7.328    0.000   57.259    0.000 cubes.py:9(all_rotations)
      8/1    8.456    1.057  108.572  108.572 cubes.py:96(generate_polycubes)
   669114    0.179    0.000    0.179    0.000 fromnumeric.py:2328(_any_dispatcher)
   669114    1.046    0.000    5.503    0.000 fromnumeric.py:2333(any)
    16302    0.004    0.000    0.004    0.000 fromnumeric.py:3241(_around_dispatcher)
     8151    0.010    0.000    0.031    0.000 fromnumeric.py:3245(around)
     8151    0.011    0.000    0.067    0.000 fromnumeric.py:3754(round_)
  2249755    1.892    0.000    4.635    0.000 fromnumeric.py:51(_wrapfunc)
  1338228    0.290    0.000    0.290    0.000 fromnumeric.py:546(_swapaxes_dispatcher)
  1338228    1.331    0.000    3.073    0.000 fromnumeric.py:550(swapaxes)
   903376    0.196    0.000    0.196    0.000 fromnumeric.py:597(_transpose_dispatcher)
   903376    1.029    0.000    3.903    0.000 fromnumeric.py:601(transpose)
   669114    1.598    0.000    4.458    0.000 fromnumeric.py:69(_wrapreduction)
   669114    0.590    0.000    0.590    0.000 fromnumeric.py:70(<dictcomp>)
  1637369    0.357    0.000    0.357    0.000 function_base.py:154(_rot90_dispatcher)
  1637369   11.148    0.000   40.823    0.000 function_base.py:158(rot90)
  1693830    0.342    0.000    0.342    0.000 function_base.py:248(_flip_dispatcher)
  1693830    8.105    0.000   16.712    0.000 function_base.py:252(flip)
  3387660    0.842    0.000    0.842    0.000 index_tricks.py:765(__getitem__)
        1    0.000    0.000    0.000    0.000 multiarray.py:1079(copyto)
  1578102    0.333    0.000    0.333    0.000 multiarray.py:1175(packbits)
  1693830    5.360    0.000    7.474    0.000 numeric.py:1348(normalize_axis_tuple)
  1693830    1.112    0.000    1.469    0.000 numeric.py:1398(<listcomp>)
        1    0.000    0.000    0.000    0.000 numeric.py:150(ones)
  1693830    0.227    0.000    0.227    0.000 {built-in method _operator.index}
     8151    0.001    0.000    0.001    0.000 {built-in method builtins.callable}
        1    0.000    0.000  108.627  108.627 {built-in method builtins.exec}
  2249755    0.394    0.000    0.394    0.000 {built-in method builtins.getattr}
  1693830    0.290    0.000    0.290    0.000 {built-in method builtins.hasattr}
  5025116    0.608    0.000    0.608    0.000 {built-in method builtins.len}
       94    0.009    0.000    0.009    0.000 {built-in method builtins.print}
       87    0.000    0.000    0.000    0.000 {built-in method builtins.round}
   903376    1.746    0.000    1.746    0.000 {built-in method numpy.arange}
   247491    0.432    0.000    0.432    0.000 {built-in method numpy.array}
  1637369    0.242    0.000    0.242    0.000 {built-in method numpy.asanyarray}
    16302    0.008    0.000    0.008    0.000 {built-in method numpy.asarray}
7844473/5230965    8.905    0.000   56.396    0.000 {built-in method numpy.core._multiarray_umath.implement_array_function}
  1693830    0.357    0.000    0.357    0.000 {built-in method numpy.core._multiarray_umath.normalize_axis_index}      
     8152    0.008    0.000    0.008    0.000 {built-in method numpy.empty}
  1355064    0.495    0.000    0.495    0.000 {method 'add' of 'set' objects}
    56461    0.012    0.000    0.012    0.000 {method 'append' of 'list' objects}
     8151    0.021    0.000    0.021    0.000 {method 'astype' of 'numpy.generic' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
  1578102    3.479    0.000    3.479    0.000 {method 'flatten' of 'numpy.ndarray' objects}
     8151    0.002    0.000    0.002    0.000 {method 'get' of 'dict' objects}
   669114    0.090    0.000    0.090    0.000 {method 'items' of 'dict' objects}
    16302    0.037    0.000    0.037    0.000 {method 'nonzero' of 'numpy.ndarray' objects}
     8151    0.008    0.000    0.008    0.000 {method 'ravel' of 'numpy.generic' objects}
     8151    0.002    0.000    0.002    0.000 {method 'ravel' of 'numpy.ndarray' objects}
   669114    2.180    0.000    2.180    0.000 {method 'reduce' of 'numpy.ufunc' objects}
     8151    0.009    0.000    0.009    0.000 {method 'round' of 'numpy.ndarray' objects}
  1338228    0.499    0.000    0.499    0.000 {method 'swapaxes' of 'numpy.ndarray' objects}
   903376    1.842    0.000    1.842    0.000 {method 'transpose' of 'numpy.ndarray' objects}

A Summary Of Problems to be solved.

pcube file format discussion

For (eventual) implementations that may support reading & writing cubes directly from disk to avoid having to store them in RAM, it would be really beneficial if we can write the count of cubes at the end of the file as well.

If writing the count to the end of the file is supported, cubes can be written to it in a streaming fashion and the count can be added at the end. If it's not supported, you have to rewrite the entire file from the beginning in order to fit the LEB128 encoding into the header.

Is there an easy way to support this?

I think we should definitely add a byte to the header that is just flags, with 1 bit explicitly reserved for increasing said header size (if we ever find more than 7 flags we need).

Polyominoes

Hello there,

I have managed to adapt the C code written by @snowmanam2 to the case of generating 2D polyominoes.
It is available in this repo.
I have also made an online viewer similar to the polycube viewer. Here is the repo.
I have computed the datasets till n=16, it was very fast like under minute or so.

I think I'll try tackling the quest of polyhypercube next, I hope my brain does not melt in the process.
I will also make the online viewer if I can find a suitable way to project the hypercubes with threejs.
Any help is appreciated.

Licence for opencubes

This was brought up for the original cubes repository and an MIT Licence was added however this repository still has no licence.
To my understanding (IANAL) we would need to get agreement from every contributor for licensing this which may have been the issue as to why it wasnt fixed at the same time but would also imply that this is something that is easier to fix sooner rather than later.

Rotation invariant representation for 2D shapes

I have been tinkering around to find a rotation invariant representation for the shapes, and so far it seems that I have found one for at least 2D. I am now trying to extend this to 3D.

My idea was inspired by Local Binary Patterns (LBP). It is used in machine vision to extract features from an image.

The idea

• See my repository for more details

To represent a shape in a rotation invariant way, the encoding should be based on relating all the cubes to all the other cubes in the shape, instead of some reference frame around the shape (because we have no way of fixing that reference frame). This makes it so that it doesn't matter in which order we compute and we still get the same result.

We first figure out a way to represent all distinct Moore neighborhoods of a cell. In 2D, there are 2^9 = 512 possible distinct Moore neighborhoods (the central cell plus all the 8-neighbors).

Then we calculate the distinct cell values for all the cells in the shape (including the empty cells) to get a feature matrix. That matrix should be unique for all shapes, but it is the same for rotated variants of that same shape.

Example in 2D

Say, we have a shape like this:

1. Add zero padding for calculation

We start by adding 2 levels of zero padding around the shape. One level is needed because the feature matrix is always larger than the shape (since all neighbors affect the cell value), and one level is needed to make sure we don't try reading outside the matrix when calculating the cell values.

2. Calculate the cell values using modified LBP

Then we can start calculating the feature matrix. We start from cell (1, 1), the second cell of the second row (so that there are not edges next to it).

To understand how the cell value is calculated, see the figure. The cell (1, 1) is highlighted in orange and its Moore neighborhood is enumerated from 1 to 9 (starting from top-left, spiraling into the center):

Since there are 512 distinct Moore neighborhoods, 9 bits is sufficient to represent all of them. We start setting the bits one by one, starting from the MSB:

1. Start from the top-left neighbor. If there is a block (=green), set the first bit (MSB) of the cell value to 1, otherwise 0.
2. Then proceed to the top-mid neighbor. Again, if there is a block, set the second bit of the cell value to 1.
3. Repeat for the top-right neighbor and set bit position 3 accordingly.
4. Repeat for the mid-right neighbor and set bit position 4 accordingly.
5. Repeat for the bottom-right neighbor and set bit position 5 accordingly.
6. Repeat for the bottom-mid neighbor and set bit position 6 accordingly.
7. Repeat for the bottom-left neighbor and set bit position 7 accordingly.
8. Repeat for the mid-left neighbor and set bit position 8 accordingly.
9. Repeat for the center cell (the cell itself) and set bit position 9 (LSB) accordingly.

Now we have a 9-bit cell value for the first cell: 00001_0000, which is 16 in decimal. We can then move to the next cell. This is repeated for all cells except cells that are on the edge (so skipping the first and last columns as well as the first and last rows; they would be zeros anyway).

Here are all the cell values for the example shape:

Then, the feature vector of this shape is:

[[ 16  24  28  12   4   0]
 [ 32  33  51  27  14   4]
 [ 64 192 496 425 263   2]
 [  0   0 112 201 390 256]
 [  0   0  96 129 258   0]
 [  0   0  64 128 256   0]]

3. Calculate the sum of the matrix values

Finally, we sum all the cell values in the feature matrix to get the fingerprint of that shape. This fingerprint is the same for all rotations of the same shape.

Comparing the variants

Here are all the variants of the example shape, along with their feature matrices and their sums (the fingerprints).

Matrix:
[[ 16  24  28  12   4   0]
 [ 32  33  51  27  14   4]
 [ 64 192 496 425 263   2]
 [  0   0 112 201 390 256]
 [  0   0  96 129 258   0]
 [  0   0  64 128 256   0]]

Sum: 3577

Matrix:
[[  0  16   8   4   0   0]
 [ 16  56  29  30  12   4]
 [ 48 105 167 291   3   2]
 [112 201 454 448 384 256]
 [ 96 129 258   0   0   0]
 [ 64 128 256   0   0   0]]

Sum: 3577

Matrix:
[[  0  16   8   4   0   0]
 [  0  48   9   6   0   0]
 [ 16 120 141 262   0   0]
 [ 32 113 155 286  12   4]
 [ 64 224 417 291   3   2]
 [  0  64 192 448 384 256]]

Sum: 3577

Matrix:
[[  0   0   0  16   8   4]
 [  0   0   0  48   9   6]
 [ 16  24  28 124 141 262]
 [ 32  33  51 107 135 258]
 [ 64 192 480 449 386 256]
 [  0   0  64 128 256   0]]

Sum: 3577

As you see, all variants have the same fingerprint 3577!

Final words

As of yet, I have not been able to prove or disprove this method, but so far it seems very promising. I have tested it on tens of various shapes in 2D, as you can see in my repository.

I want to prove that this method actually works and never gives the same fingerprint to 2 different shapes, then proceed to applying this to 3D, if possible.

Javascript Online Viewer

Hello there,
I have made an application so that you may load and visualize a polycube dataset in your browser.

You can load dataset from your computer in a text file format where each line is xyz111000110001...
I have seen the post about the pcube format but it was too complex for me to implement fast in Javascript.
So I also propose to download polycubes from 3 to 9 in the correct format as an archive (if you trust it):
or individually polycubes 3, polycubes 4, polycubes 5, polycubes 6, polycubes 7, polycubes 8, polycubes 9.

Controls:

• Mouse lets you zoom in, zoom out and rotate around.
• Up/Down arrows are for moving to the previous/next polycube
• Left/Right arrows are for moving to theprevious/next polycube by the step indicated by the slider (default is 10^3=1000)
  You can change the step value using the slider and it is in log10 scale so that you can jump to any index in a large dataset quickly.
• Key 5 lets you reset the camera position

I also spend some part of my weekend working on the generation of all polycubes that fits across a 3x3x3 cube:
there are 1 551 811 of them, took 20min on my computer to compute with some Python script that I have optimized a bit since it was 90min at first.
Here are the dataset: zip archive (3.3MB) and txt file (50MB).

I hope you will find it useful.
There are many improvements that can be brought to this viewer, let me know if you think of good ones.

Edit 03/09/2023 : All linked files have been removed.

Solution where memory usage isn't an issue

@hidny has come up with a solution where memory usage isn't the limitation. It also appears faster. You can find my implementation in https://gitlab.com/dzamlo/polycubes2

Unfortunaly, I cannot make a pull request because my laptop died.