Cubester - Rubik's Cube Simulator

During 2022, my son has played with Rubik's cubes, just as I did in his age. He has solved the original 3 by 3 cube in 15 seconds (while I never solved it - although back then there were no YouTube instruction videos). He can also solve the 2 by 2, 4 by 4, 5 by 5 and the Megaminx. So I figured it would be fun to model a cube in Python as well. And yes, my son can solve it using the program as well. But perhaps the more physical user interface of having the actual plastic cube is better - it took him about 10 minutes for a 3 by 3 cube...

This implementation can model cubes of almost any reasonable size. As usual, I used Python, NumPy, and pygame to do this, and all code can be found in my GitHub repository.

You can set the cube to Auto Rotate or be still; rotate it around any of its axes; use automatic shuffling or undo; and of course turn the individual discs (layers) of cubies (small cubes) as in a real cube. The buttons on the left can be operated by using the mouse or the shortcut keys.

Cubic Rubik Modelling

There are many modelling decisions involved when modelling a Rubik's Cube. I wanted to achieve the following:

• model any size (i.e. any number of cubies (small cubes) per cube side)
• enable full 3D rotation
• enable animated rotation of any disc (or layer) of cubies

So I created a class Cube to hold the data. The main parameter when adding a new cube is size i.e. the number of cubies per side, being 3 for the original cube.

        self.cube = -np.ones((size, size, size), dtype=np.int32)
        self.cubies_nr = size ** 3 - max(0, size - 2) ** 3  # cubies_nr = outer edge small cubes; subtract inner cubes from all cubes.
        self.nodes = np.array([
            [-1, -1, -1],
            [-1, -1, +1],
            [+1, -1, +1],
            [+1, -1, -1],
            [-1, +1, -1],
            [-1, +1, +1],
            [+1, +1, +1],
            [+1, +1, -1]
            ])
        self.surfaces = np.array([
            [0, 1, 2, 3],
            [7, 6, 5, 4],
            [0, 3, 7, 4],
            [1, 0, 4, 5],
            [2, 1, 5, 6],
            [3, 2, 6, 7]
            ])
        # colors of cube sides (1-6) and "non-colored" side (color 0)
        self.color_0 = 60
        self.colors = np.array([
            [self.color_0, self.color_0, self.color_0],
            [255, 255, 255],
            [255, 255,   0],
            [  0, 200,   0],
            [  0,   0, 200],
            [235, 150,   0],
            [200,   0,   0]
            ], dtype=np.uint8)

First I set up the cube as having three dimensions of size size, each item holding a reference -1 to cubies. The number of cubies required is not the full size^3, as the cubies below the surface do not really exists - as in physical cubes. In a 3 by 3 cube, there are 3^3 = 27 cubies less the 1 inside = 26 cubies, while for e.g. in a 6 by 6 cube there are 6^3 - 4^3 = 152 cubies (instead of 216). Then I define the "standard" nodes (the eight corners of each cube), the surfaces (six flat surfaces with four corners each), and the colors (six cube colors and the cubie "non-colored" side color, in RGB).

Not shown here, I then go through each cubie and define its 3D coordinates with respect to the Cube center, and the color of each surface. 

When building the class, I originally designed it thinking about a 3 by 3 cube - a very small one, that is. It turns out the class is not super efficient for big cubes (like greater than 30 by 30), for which a different structure for data might be better suited. However, the huge cubes are anyway more of a curiosity, perhaps interesting conceptually, but rather difficult to work with.

Rotating the Rotated

Now with the 3D coordinates and surfaces defined for each cubie, we can draw it on screen. First it can of course be rotated in 3D. This I have covered in my previous posts, so will not do that again - instead, I will go through how to rotate this 3D object around any selected vector in 3D space. While the rotation angles self.angles in def rotate are defined in X, Y, and Z rotation, they have a "priority order": first X, then Y, then Z (this could be any order, but I have chosen this one). This means that if the angles are all zero, the object can be rotated around any of the three axis by altering the respective angle. But, if the angles are something else (the object has already been rotated), changing the angle on one axis does not rotate the object around the original axis, except for the X-axis, which comes first. You can think of it as follows: First we rotate around the X axis. Then, the already rotated object is rotated around the Y axis, and the resulting rotated object then again around the Z axis.

Rotating around a cube axis means that when we have already rotated the cube to some random angle (x, y, z), we wish to rotate it around one of the original axes. This is the code:

    def angle_add(self, cube, time, prev_time):

        unit_vector = np.array([0, 0, 0], dtype=np.float)
        time_adjustment = (time - prev_time) / (20 * 17)
        keys = pygame.key.get_pressed()
        if keys[pygame.K_RIGHT] or self.mouse_rotate == 'RIGHT':
            unit_vector[0] += 1
        if keys[pygame.K_LEFT] or self.mouse_rotate == 'LEFT':
            unit_vector[0] -= 1
        if keys[pygame.K_UP] or self.mouse_rotate == 'UP':
            unit_vector[1] += 1
        if keys[pygame.K_DOWN] or self.mouse_rotate == 'DOWN':
            unit_vector[1] -= 1
        if keys[pygame.K_PAGEUP] or self.mouse_rotate == 'PAGE UP':
            unit_vector[2] += 1
        if keys[pygame.K_PAGEDOWN] or self.mouse_rotate == 'PAGE DOWN':
            unit_vector[2] -= 1
        if not np.all(np.equal(unit_vector, np.array([0, 0, 0]))):

            # rotate along one of the cube's axes, based on the unit vector chosen. See https://en.wikipedia.org/wiki/Rotation_matrix
            # if more than one axes chosen (e.g. UP and LEFT), that will work as well.
            # first rotate the unit vector to point according to current cube angles, then rotate around it.
            unit_vector /= np.sqrt(np.sum(unit_vector ** 2))  # scale to unit vector if multiple axes used
            matrix = self.rotate_matrix(cube.angles)
            (x, y, z) = np.matmul(unit_vector, matrix)
            angle = 1.0 * time_adjustment
            sa = np.sin(angle)
            ca = np.cos(angle)
            matrix_uv = np.array([[ca + x * y * (1 - ca)    , x * y * (1 - ca) - z * sa, x * z * (1 - ca) + y * sa],
                                  [y * x * (1 - ca) + z * sa, ca + y * y * (1 - ca)    , y * z * (1 - ca) - x * sa],
                                  [z * x * (1 - ca) - y * sa, z * y * (1 - ca) + x * sa, ca + z * z * (1 - ca)    ]])

            # final rotation matrix is current matrix rotated around unit vector.
            matrix_final = np.matmul(matrix, matrix_uv)
            # and the resulting angles can be calculated from that matrix.
            angle_x = np.arctan2(-matrix_final[1, 2], matrix_final[2, 2])
            angle_y = np.arctan2(matrix_final[0, 2], np.sqrt(1 - matrix_final[0, 2] ** 2))
            angle_z = np.arctan2(-matrix_final[0, 1], matrix_final[0, 0])
            cube.angles = np.array([angle_x, angle_y, angle_z])

There are many ways to achieve the desired effect, but this the most general, being able to rotate the object around any vector in 3D space going through the origo, and add to it a general rotation around the X, Y, Z as described.

First I get the rotation from the user; specific keys are checked and the respective rotation added to the unit_vector. The self.mouse_rotate is set elsewhere if the mouse pointer is on top of one of the rotation buttons displayed, and the left mouse button is pressed - rotation can be activated by either the keys or the mouse. Note that all the axes are additive, so if both left and up cursor keys are pressed, the resulting vector is [-1, 1, 0]. This is also why a scaling (dividing the vector by its length) is needed to actually make this vector a unit vector (a vector with a length of 1): unit_vector /= np.sqrt(np.sum(unit_vector ** 2)).

Then, I need to get the normal rotation matrix, to rotate the unit vector with it. Now (x, y, z) represent the unit vector rotated with the cube angles. Then, selecting some angle = 1.0 * time_adjustment representing the speed of rotation around that unit vector, I can calculate another rotation matrix matrix_uv using the rotated unit vector, and rotate it with the original matrix to produce a combined, final rotation matrix: matrix_final = np.matmul(matrix, matrix_uv). This would actually be enough to rotate the cube, but I can also go "backwards" a bit and get the respective angles for that rotation in X, Y, and Z (again the order is important), and store them as new cube.angles. This helps to keep track of where the current rotation is at any given moment.

Solving and Shuffling

To be interesting, the cube needs to work. It must be possible to rotate a single disc (or layer) as in a physical cube. That is why, in defining the cube, I have given all cubies their own, independent nodes. This way, a group of cubies can be rotated 90 degrees to either direction to emulate turning a disc (or layer) on a physical cube. There are two ways to handle what happens after such a turn in the structure I have chosen: 1) the cubies can remain where they are, after the turn, and the data in the cube, representing the cubie in each of the cube positions, can be updated accordingly; or 2) the cubies could be returned to their original positions, but their properties (surface colors) could be changed so that the cube looks as if they were moved. I chose the (1).

The animated turning of any single disc is realized by rotating the required cubies first around the required axis, and then using the resulting rotated coordinates as a basis when rotating the whole cube to its desired angles.

For a nice user experience, there is an automatic shuffling option, randomly turning the discs. And to help aspiring cube masters, the keys for turning any specific disc are shown next to the corresponding disc. I have used nine keys for each of the three dimensions:

        self.disc_keys_x = 'qwertyuio'
        self.disc_keys_y = 'asdfghjkl'
        self.disc_keys_z = 'zxcvbnm,.'

You should note that these follow my Finnish keyboard and should be changed for any other keyboard layout. I was unable to figure out a way to get key positions instead of key codes. Of course, for a 3 by 3 cube, only the first three keys for each axes can be used, and similarly for other sizes. When cube size grows beyond the key list length (here 9 keys), CTRL and ALT keys can be used to rotate higher level discs. The SHIFT key can be used to reverse the direction of any turn, until the cube size grows so big that it, too, is needed to get to higher levels. (Note that at least Windows can "steal" some key combinations, if they are defined as shortcuts; my HP computer had CTRL+ALT+S for some HP help program, which needed to be removed.)