A MILP based Happy Cube©® Solver

Posted on 31. August 202415. September 2024 by Fabian Leuthold

In this post, we’re developing a Happy Cube solver using a mixed integer linear program (MILP). Previously, we tackled this problem on a micro:bit using dynamic programming and backtracking. Now, we’re trying a different approach by creating a mathematical model to solve this combinatorial challenge.

The idea is to structure the MILP so that it tells us which piece goes where in the cube’s net, how many 90-degree counterclockwise rotations are needed, and whether the piece needs to be flipped horizontally before placement to swap the front and back sides.

Our first Test Instance

Before we start modeling, we’ll choose a test instance to validate our model. We’ll lay out the pieces of a cube in front of us and number them for reference.

To read the information, I’ve come up with a simple file format. Similar to the micro:bit program, we’ll save the binary pattern of the edges for each piece sequentially. We start by listing the top edge first, then move clockwise around the edges. So, the first 5 bits represent the top edge, the next 5 bits are for the right edge, followed by 5 bits for the bottom edge, and finally, 5 bits for the left edge.

For clarity, I’ve decided to list the corners twice and accept some redundancy to make the model easier to interpret. The bit sequences for each line in the file are separated by semicolons to help visually organize the data and prevent errors during manual entry. These semicolons are ignored when reading the file. 😜

The following text file – also available here – shows the representation of the problem as depicted in the photo above:

NAME: HC_BQ_001
TYPE: HC_BQ_INSTANCE
COMMENT: Build 3D-Cube out of 6 Happy Cube (R) tiles
NUM_TILES: 6
TILE_ID_BORDER_PATTERN
1 [0, 1, 0, 1, 0; 0, 0, 1, 0, 1; 1, 1, 0, 1, 1; 1, 1, 0, 1, 0]
2 [1, 1, 0, 1, 1; 1, 0, 1, 0, 0; 0, 0, 1, 0, 1; 1, 1, 0, 1, 1]
3 [1, 0, 1, 0, 0; 0, 0, 1, 0, 0; 0, 0, 1, 0, 1; 1, 1, 0, 1, 1]
4 [0, 1, 0, 1, 0; 0, 0, 1, 0, 0; 0, 0, 1, 0, 0; 0, 1, 0, 1, 0]
5 [0, 1, 0, 1, 0; 0, 0, 1, 0, 0; 0, 1, 0, 1, 1; 1, 0, 1, 0, 0]
6 [0, 1, 0, 1, 0; 0, 0, 1, 0, 0; 0, 1, 0, 1, 0; 0, 0, 1, 0, 0]
EOF

The following image shows how the bit pattern for the row of Tile-ID 1 in the text file is formed:

*Visualization of the generation of the pit pattern for tile with ID (1) based on the tile’s position & orientation*

Some Ideas for the Model

The model might use a preprocessed parameter to know the border of the tiles of all individual pieces in all possible rotation- and swap-states. The parameter can be determined in preprocessing and later multiplied with the decision variable, which virtually allows the model to test all possible configurations.
This allows it to determine if two pieces of a given rotation/flip configuration fit together or not. Wherever one piece sticks out, there needs to be a gap in the counterpart piece.

*three bits center bits determine, whether two parts fit together*

When you look closely at the cubes, you’ll notice that the corners are special. In each corner, three pieces meet, but only one of them is allowed to fill the corner.

*every corner results from the sum of three single bits: here empty corner (disallowed!)*

The following graphic provides an overview of the corners and edges in the cube net and shows which pieces make up each corner and which ones combine to the same edge.

*The blue colored dots mark identic corners, the orange ones mark identic edges in the cube net*

Modeling the Problem as a MILP

We’ll begin by modeling the MILP following a standard approach and see how far we get. The process includes these steps:

Define decision variables
Define auxiliary variables
Introduce indices / sets
Identify relevant parameters
Derive constraints
Derive the objective function
Validate the model

Define decision Variables

To determine the structure of the decision variables, it helps to clarify what we want to extract from the solution.

As mentioned earlier, we expect the model to specify which piece goes into each position of the cube’s net. Additionally, we need to know if the piece should be flipped horizontally before placement and whether it needs to be rotated 0, 1, 2, or 3 times by 90 degrees counterclockwise before fitting it into the correct spot.

It probably makes sense to first think about the indices and sets, as this will help us better define the decision variables. So, we’ll jump to step 3 of our MILP recipe and then circle back to step 1.

Introduce Indices / Sets

First, we need an ID to name the parts on the table. We use the letter $i$ for this. Each index $i$ represents one of the six puzzle pieces, $i \in I =\{1, 2, 3, 4, 5, 6\}$ .
We need an ID to tell, where the pieces must go in the cube net, in order to know, how we should put the pieces together to form a cube. Let’s take letter $j$ for this. Each index $j$ represents one of the possible positions, a piece could go in the cube’s net: $j \in J =\{1, 2, 3, 4, 5, 6\}$ .
We need another ID to tell us the orientation, 1 meaning, we have to flip a piece horizontally, before placing it in the cube net, 0 meaning we don’t have to flip it. We introduce another index and use $k \in K = \{0, 1\}$ .
Finally, we need an index, telling, how many times a given piece must be rotated counter clockwise before being placed in the cube net. We introduce a last index and use letter $l \in L=\{0, 1, 2, 3\}$ .

*Position indices of parts in the cube net, represented by index j*

Define decision Variables, the 2nd

Let’s get back to our question: what should we name our decision variable? We could use this:

$X_{ijkl} = 1$ , exactly if piece $i$ needs to be flipped horizontally if $k=1$ and rotated counterclockwise $l$ times in 90 degree steps before being placed in position $j$ of the cube net.

Since we don’t need any auxiliary variables right now, let’s move straight to step 4.

Identify relevant Parameters

Now we need a parameter to represent the binary patterns of the edges of pieces based on their position and orientation. We’ll use the letter p (for “part”).

$P_{iklm} = 1$ , if the border pattern of part $i$ at index $m \in M={1, 2, ..., 20}$ is 1 and the part is being flipped horizontally ( $k= 1$ ) or not flipped ( $k=0$ ) and rotated counter clockwise $l$ times.

We just realized, that we need another index to choose the position of a specific bit in the bit-string of each part. For this we introduced index $m \in M={1, 2, ..., 20}$ .

Now we should be able to build our model.

Derive constraints

First we need some constraints to prevent the solver from choosing invalid variable configurations. We must make sure, that every part ends up being placed at exactly one possible position and orientation:

C1: Any piece must exactly be placed once
$\forall i \in I: \sum_{(j,k,l)} X_{i,j,k,l} = 1$
Make sure that every part goes exactly to one position in a specific orientation and rotation.

Now, let’s formulate the constraints to ensure that the parts placed in specific locations on the cube net fit together correctly, ultimately allowing them to be assembled into a 3D cube.

C2: Edge 1 (Intersection Part 1 & Part 2)
$\forall m \in \{12,13,14\}: \sum_{(i,k,l)} X_{i,1,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,2,k,l} \cdot P_{i,k,l,16 - m}$
Make sure that the bottom side of the part chosen for (cube net) location 1 fits together with the top side of the part chosen for (cube net) location 2.
C3: Edge 2 (Intersection Part 1 & Part 4)
$\forall m \in \{7,8,9\}: \sum_{(i,k,l)} X_{i,1,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,4,k,l} \cdot P_{i,k,l,11-m}$
Make sure that the right side of the part chosen for location 1 fits together with the top side of the part chosen for location 4.
C4: Edge 3 (Intersection Part 1 & Part 3)
$\forall m \in \{17,18,19\}: \sum_{(i,k,l)} X_{i,1,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,3,k,l} \cdot P_{i,k,l,21-m}$
Make sure that the left side of the part chosen for location 1 fits together with the top side of the part chosen for location 3.
C5: Edge 4 (Intersection Part 1 & Part 6)
$\forall m \in \{2, 3, 4\}: \sum_{(i,k,l)} X_{i,1,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,6,k,l} \cdot P_{i,k,l,16-m}$
Make sure that the top side of the part chosen for location 1 fits together with the bottom side of the part chosen for location 6.
C6: Edge 5 (Intersection Part 5 & Part 2)
$\forall m \in \{2, 3, 4\}: \sum_{(i,k,l)} X_{i,5,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,2,k,l} \cdot P_{i,k,l,16-m}$
Make sure that the top side of the part chosen for location 5 fits together with the bottom side of the part chosen for location 2.
C7: Edge 6(Intersection Part 5 & Part 4)
$\forall m \in \{7,8,9\}: \sum_{(i,k,l)} X_{i,5,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,4,k,l} \cdot P_{i,k,l,21-m}$
Make sure that the right side of the part chosen for location 5 fits together with the bottom side of the part chosen for location 4.
C8: Edge 7 (Intersection Part 5 & Part 6)
$\forall m \in \{12, 13, 14\}: \sum_{(i,k,l)} X_{i,5,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,6,k,l} \cdot P_{i,k,l,16-m}$
Make sure that the bottom side of the part chosen for location 5 fits together with the top side of the part chosen for location 6.
C9: Edge 8 (Intersection Part 5 & Part 3)
$\forall m \in \{17,18,19\}: \sum_{(i,k,l)} X_{i,5,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,3,k,l} \cdot P_{i,k,l,31-m}$
Make sure that the left side of the part chosen for location 5 fits together with the bottom side of the part chosen for location 3.
C10: Edge 9 (Intersection Part 6 & Part 3)
$\forall m \in \{17,18,19\}: \sum_{(i,k,l)} X_{i,6,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,3,k,l} \cdot P_{i,k,l,36-m}$
Make sure that the left side of the part chosen for location 6 fits together with the left side of the part chosen for location 3.
C11: Edge 10 (Intersection Part 6 & Part 4)
$\forall m \in \{7,8,9\}: \sum_{(i,k,l)} X_{i,6,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,4,k,l} \cdot P_{i,k,l,16-m}$
Make sure that the right side of the part chosen for location 6 fits together with the right side of the part chosen for location 4.
C12: Edge 11 (Intersection Part 2 & Part 3)
$\forall m \in \{17,18,19\}: \sum_{(i,k,l)} X_{i,2,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,3,k,l} \cdot P_{i,k,l,26-m}$
Make sure that the left side of the part chosen for location 2 fits together with the right side of the part chosen for location 3.
C13: Edge 12 (Intersection Part 2 & Part 4)
$\forall m \in \{7,8,9\}: \sum_{(i,k,l)} X_{i,2,k,l} \cdot P_{i,k,l,m} = 1 - \sum_{(i,k,l)} X_{i,4,k,l} \cdot P_{i,k,l,26-m}$
Make sure that the right side of the part chosen for location 2 fits together with the left side of the part chosen for location 4.

When examining the cube’s corners in relation to our model, we noticed that the corners are a special case since three pieces meet there. The model needs to ensure that only one of the three pieces fills the corner. We’re now trying to ensure this with the following constraints for the 8 corners:

C14: Corner 1 (Intersection of Parts 1, 2, 3)
$\sum_{(i,k,l)} X_{i,1,k,l} \cdot P_{i,k,l,15} + \sum_{(i,k,l)} X_{i,2,k,l} \cdot P_{i,k,l,20} + \sum_{(i,k,l)} X_{i,3,k,l} \cdot P_{i,k,l,5} = 1$
Make sure that in corner 1 only one part of the parts chosen for location 1,2 & 3 provides a bit to fill the gap.
C15: Corner 2 (Intersection of Parts 1, 2, 4)
$\sum_{(i,k,l)} X_{i,1,k,l} \cdot P_{i,k,l,10} + \sum_{(i,k,l)} X_{i,2,k,l} \cdot P_{i,k,l,5} + \sum_{(i,k,l)} X_{i,4,k,l} \cdot P_{i,k,l,1} = 1$
Make sure that in corner 2 only one part of the parts chosen for location 1,2 & 4 provides a bit to fill the gap.
C16: Corner 3 (Intersection of Parts 2, 4, 5)
$\sum_{(i,k,l)} X_{i,2,k,l} \cdot P_{i,k,l,10} + \sum_{(i,k,l)} X_{i,4,k,l} \cdot P_{i,k,l,15} + \sum_{(i,k,l)} X_{i,5,k,l} \cdot P_{i,k,l,5} = 1$
Make sure that in corner 3 only one part of the parts chosen for location 2, 4 & 5 provides a bit to fill the gap.
C17: Corner 4 (Intersection of Parts 2, 3, 5)
$\sum_{(i,k,l)} X_{i,2,k,l} \cdot P_{i,k,l,15} + \sum_{(i,k,l)} X_{i,3,k,l} \cdot P_{i,k,l,10} + \sum_{(i,k,l)} X_{i,5,k,l} \cdot P_{i,k,l,1} = 1$
Make sure that in corner 4 only one part of the parts chosen for location 2, 3 & 5 provides a bit to fill the gap.
C18: Corner 5 (Intersection of Parts 3, 5, 6)
$\sum_{(i,k,l)} X_{i,3,k,l} \cdot P_{i,k,l,15} + \sum_{(i,k,l)} X_{i,5,k,l} \cdot P_{i,k,l,15} + \sum_{(i,k,l)} X_{i,6,k,l} \cdot P_{i,k,l,1} = 1$
Make sure that in corner 5 only one part of the parts chosen for location 3, 5 & 6 provides a bit to fill the gap.
C19: Corner 6 (Intersection of Parts 4, 5, 6)
$\sum_{(i,k,l)} X_{i,4,k,l} \cdot P_{i,k,l,10} + \sum_{(i,k,l)} X_{i,5,k,l} \cdot P_{i,k,l,10} + \sum_{(i,k,l)} X_{i,6,k,l} \cdot P_{i,k,l,5} = 1$
Make sure that in corner 6 only one part of the parts chosen for location 4, 5 & 6 provides a bit to fill the gap.
C20: Corner 7 (Intersection of Parts 1, 4, 6)
$\sum_{(i,k,l)} X_{i,1,k,l} \cdot P_{i,k,l,5} + \sum_{(i,k,l)} X_{i,4,k,l} \cdot P_{i,k,l,5} + \sum_{(i,k,l)} X_{i,6,k,l} \cdot P_{i,k,l,10} = 1$
Make sure that in corner 7 only one part of the parts chosen for location 1, 4 & 6 provides a bit to fill the gap.
C21: Corner 8 (Intersection of Parts 1, 3, 6)
$\sum_{(i,k,l)} X_{i,1,k,l} \cdot P_{i,k,l,1} + \sum_{(i,k,l)} X_{i,3,k,l} \cdot P_{i,k,l,1} + \sum_{(i,k,l)} X_{i,6,k,l} \cdot P_{i,k,l,15} = 1$
Make sure that in corner 8 only one part of the parts chosen for location 1, 3 & 6 provides a bit to fill the gap.

Derive the Objective Function

This problem is essentially a feasibility task, so I initially set the objective function to $min \; 0$ .

However, while solving different cubes, I noticed the model sometimes performs unnecessary flips with symmetrical pieces or rotates them by 180 degrees, even though this doesn’t change the situation.

To address this, I created an auxiliary variable to count the number of flips and 180-degree rotations and then minimized it. So, we need an additional 22nd constraint:

C22: Memorize the number of horizontal flips and 2-turns
$Y \ge \sum_{(i,j,l)} X_{i,j,1,l} + \sum_{(i,j,k)} X_{i,j,k,2}$
Collect the sum of the number of horizontal flips and the number of 180° rotations in Y. Since we minimize the objective function containing $Y$ , the solver would never choose this number to exceed the actual sum of these two values.

Now we can minimize the auxiliary variable, which in turn minimizes the number of swaps and 180-degree rotations:

$min \;Y$

The Implementation

Here you can find the complete source:

# Happy Cube (C)(R) Solver based on MILP
# Author: Fabian Leuthold
# Date: 2024 - 08 - 31

import pytest
import os
import numpy as np
from ortools.linear_solver import pywraplp


def read_hcq_file(file):
    """
    reads a happy cube (C)(R) instance from a text-based hcq-file

    :param file: the file path
    :return: a list of lists containing the binary border-pattern of every tile as they lay horizontally placed and
             enumerated by integers from 1 to 6.
    """
    with open(file, 'r') as f:
        lines = f.readlines()

    binary_lists = []

    for line in lines:
        # Check if line contains a tile pattern
        if '[' in line and ']' in line:
            # Extract the part within the brackets
            start = line.index('[') + 1
            end = line.index(']')
            pattern_str = line[start:end]
            # Replace all semicolons with commas
            pattern_str = pattern_str.replace(';', ',')
            # Split by commas, then convert to integers
            pattern = [int(num) for num in pattern_str.split(',')]
            assert len(pattern) == 5 * 4, "Expected each partern to have 20 binary values but failed"
            binary_lists.append(pattern)

    assert len(binary_lists) == 6, "Expected to find 1 border-pattern for each cube tile but failed"

    return binary_lists


def convert_4d_var(var, I, J, K, L):
    """
    converts the solver X in a np-X

    :param var: X_ijkl
    :param I: list of i
    :param J: list of j
    :param K: list of k
    :param L: list of l
    :return: np X
    """
    retval = dict()
    for i in I:
        for j in J:
            for k in K:
                for l in L:
                    if var[i, j, k, l].solution_value():
                        str_swap = ", swap it horizontally," if k else ""
                        str_rot = f" rotate it {l} times 90 degrees counter-clockwise" if l != 0 else ""
                        retval[j] = f"Take part {i + 1}{str_swap}{str_rot} and put it on place {j + 1}."
    return retval


def solve_hcq_milp(hcq_tiles):
    # Init return values
    solution = []

    # Create the model.
    solver_name = 'CP-SAT'
    # solver_name = 'SCIP'
    # solver_name = 'GUROBI_MIP'
    solver = pywraplp.Solver.CreateSolver(solver_name)
    # solver.EnableOutput()

    # define index sets
    I = range(len(hcq_tiles))  # the ID (1-6) of the tiles as they are currently placed on the table
    J = range(len(hcq_tiles))  # the position ID (1-6) as they shall be placed within the cube pattern;
    # 1, 2, 5, 6 (vertically), 3, 2, 4 (horizontally)
    K = range(2)  # apply horizontal swap to tile before placing in the net (=1) or not (=0)?
    L = range(4)  # rotate tile counterclockwise for 0, 1, 2, 3 x 90° ?
    M = range(20)  # bit-string length for every part (4*5)

    # define parameters
    P = do_preprocessing(hcq_tiles, I, K, L, M)

    # declare decision variables
    X = {}
    for i in I:
        for j in J:
            for k in K:
                for l in L:
                    X[i, j, k, l] = solver.BoolVar('item_i%ij%ik%il%i' % (i, j, k, l))

    Y = solver.IntVar(0, 6, "Y")

    # C1: Any piece must exactly be placed once
    for i in I:
        solver.Add(sum(X[i, j, k, l] for j in J for k in K for l in L) == 1)
    # C2: Edge 1 (Intersection Part 1 & Part 2)
    for m in [11, 12, 13]:
        solver.Add(sum(X[i, 0, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 1, k, l] * P[i, k, l, 14 - m] for i in I for k in K for l in L))
    # C3: Edge 2 (Intersection Part 1 & Part 4)
    for m in [6, 7, 8]:
        solver.Add(sum(X[i, 0, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 3, k, l] * P[i, k, l, 9 - m] for i in I for k in K for l in L))
    # C4: Edge 3 (Intersection Part 1 & Part 3)
    for m in [16, 17, 18]:
        solver.Add(sum(X[i, 0, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 2, k, l] * P[i, k, l, 19 - m] for i in I for k in K for l in L))
    # C5: Edge 4 (Intersection Part 1 & Part 6)
    for m in [1, 2, 3]:
        solver.Add(sum(X[i, 0, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 5, k, l] * P[i, k, l, 14 - m] for i in I for k in K for l in L))
    # C6: Edge 5 (Intersection Part 5 & Part 2)
    for m in [1, 2, 3]:
        solver.Add(sum(X[i, 4, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 1, k, l] * P[i, k, l, 14 - m] for i in I for k in K for l in L))
    # C7: Edge 6(Intersection Part 5 & Part 4)
    for m in [6, 7, 8]:
        solver.Add(sum(X[i, 4, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 3, k, l] * P[i, k, l, 19 - m] for i in I for k in K for l in L))
    # C8: Edge 7 (Intersection Part 5 & Part 6)
    for m in [11, 12, 13]:
        solver.Add(sum(X[i, 4, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 5, k, l] * P[i, k, l, 14 - m] for i in I for k in K for l in L))
    # C9: Edge 8 (Intersection Part 5 & Part 3)
    for m in [16, 17, 18]:
        solver.Add(sum(X[i, 4, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 2, k, l] * P[i, k, l, 29 - m] for i in I for k in K for l in L))
    # C10: Edge 9 (Intersection Part 6 & Part 3)
    for m in [16, 17, 18]:
        solver.Add(sum(X[i, 5, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 2, k, l] * P[i, k, l, 34 - m] for i in I for k in K for l in L))
    # C11: Edge 10 (Intersection Part 6 & Part 4)
    for m in [6, 7, 8]:
        solver.Add(sum(X[i, 5, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 3, k, l] * P[i, k, l, 14 - m] for i in I for k in K for l in L))
    # C12: Edge 11 (Intersection Part 2 & Part 3)
    for m in [16, 17, 18]:
        solver.Add(sum(X[i, 1, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 2, k, l] * P[i, k, l, 24 - m] for i in I for k in K for l in L))
    # C13: Edge 12 (Intersection Part 2 & Part 4)
    for m in [6, 7, 8]:
        solver.Add(sum(X[i, 1, k, l] * P[i, k, l, m] for i in I for k in K for l in L) == 1 - sum(
            X[i, 3, k, l] * P[i, k, l, 24 - m] for i in I for k in K for l in L))
    # C14: Corner 1 (Intersection of Parts 1, 2, 3)
    solver.Add(sum(X[i, 0, k, l] * P[i, k, l, 14] for i in I for k in K for l in L) +
               sum(X[i, 1, k, l] * P[i, k, l, 19] for i in I for k in K for l in L) +
               sum(X[i, 2, k, l] * P[i, k, l, 4] for i in I for k in K for l in L) == 1)
    # C15: Corner 2 (Intersection of Parts 1, 2, 4)
    solver.Add(sum(X[i, 0, k, l] * P[i, k, l, 9] for i in I for k in K for l in L) +
               sum(X[i, 1, k, l] * P[i, k, l, 4] for i in I for k in K for l in L) +
               sum(X[i, 3, k, l] * P[i, k, l, 0] for i in I for k in K for l in L) == 1)
    # C16: Corner 3 (Intersection of Parts 2, 4, 5)
    solver.Add(sum(X[i, 1, k, l] * P[i, k, l, 9] for i in I for k in K for l in L) +
               sum(X[i, 3, k, l] * P[i, k, l, 14] for i in I for k in K for l in L) +
               sum(X[i, 4, k, l] * P[i, k, l, 4] for i in I for k in K for l in L) == 1)
    # C17: Corner 4 (Intersection of Parts 2, 3, 5)
    solver.Add(sum(X[i, 1, k, l] * P[i, k, l, 14] for i in I for k in K for l in L) +
               sum(X[i, 2, k, l] * P[i, k, l, 9] for i in I for k in K for l in L) +
               sum(X[i, 4, k, l] * P[i, k, l, 0] for i in I for k in K for l in L) == 1)
    # C18: Corner 5 (Intersection of Parts 3, 5, 6)
    solver.Add(sum(X[i, 2, k, l] * P[i, k, l, 14] for i in I for k in K for l in L) +
               sum(X[i, 4, k, l] * P[i, k, l, 14] for i in I for k in K for l in L) +
               sum(X[i, 5, k, l] * P[i, k, l, 0] for i in I for k in K for l in L) == 1)
    # C19: Corner 6 (Intersection of Parts 4, 5, 6)
    solver.Add(sum(X[i, 3, k, l] * P[i, k, l, 9] for i in I for k in K for l in L) +
               sum(X[i, 4, k, l] * P[i, k, l, 9] for i in I for k in K for l in L) +
               sum(X[i, 5, k, l] * P[i, k, l, 4] for i in I for k in K for l in L) == 1)
    # C20: Corner 7 (Intersection of Parts 1, 4, 6)
    solver.Add(sum(X[i, 0, k, l] * P[i, k, l, 4] for i in I for k in K for l in L) +
               sum(X[i, 3, k, l] * P[i, k, l, 4] for i in I for k in K for l in L) +
               sum(X[i, 5, k, l] * P[i, k, l, 9] for i in I for k in K for l in L) == 1)
    # C21: Corner 8 (Intersection of Parts 1, 3, 6)
    solver.Add(sum(X[i, 0, k, l] * P[i, k, l, 0] for i in I for k in K for l in L) +
               sum(X[i, 2, k, l] * P[i, k, l, 0] for i in I for k in K for l in L) +
               sum(X[i, 5, k, l] * P[i, k, l, 14] for i in I for k in K for l in L) == 1)

    # C22: Memorize no of flips in Y and rotations of 2
    solver.Add(
        Y >= sum(X[i, j, 1, l] for i in I for j in J for l in L) + sum(X[i, j, k, 2] for i in I for j in J for k in K))

    # objective function
    solver.Minimize(Y)

    # solve the model
    status = solver.Solve()

    if status == pywraplp.Solver.OPTIMAL:
        print("Optimal solution found.")
        solution = convert_4d_var(X, I, J, K, L)
    else:
        print("The problem does not have an optimal solution.")
    return solution


def apply_swap_operation(cur_bit_stream):
    """
    takes a bit stream from a part (20 bits in length expected) and applies the transformation which would
    occurr if the part was flipped horizontally.

    :param cur_bit_stream: the current bit stream of a given tile
    :return: the transformed bit stream (horiz. flip)
    """
    return cur_bit_stream[4::-1] + cur_bit_stream[20:14:-1] + cur_bit_stream[14:9:-1] + cur_bit_stream[9:4:-1]


def apply_rotation(cur_bit_stream, l):
    """
    rotate the part contained in cur_cur_bit_stream by applying l-times a 90 degree counter clockwise rotation

    :param cur_bit_stream: the bit stream of a given tile
    :param l: the number of counter clockwise rotations in steps of 90 degrees
    :return: the transformed bit stream
    """
    # clone list to prevent manipulating the origin data structure
    for _ in range(1, l + 1):
        cur_bit_stream = cur_bit_stream[5:10] + cur_bit_stream[10:15] + cur_bit_stream[15:20] + cur_bit_stream[:5]
    return cur_bit_stream


def get_binary_value_from_stream(hcq_tiles, i, k, l, m):
    """
    determines the binary value for a specific tile i at a given position m, swapped k or rotated l.

    :param hcq_tiles: the bit stream for all tiles as a list of lists
    :param i: the index of the tile
    :param k: the swap state (0 = not swapped, 1 = swapped)
    :param l: the rotation state (0 = not rotated, 1,2,3 = rot. in 90 degs 1,2,3 times in counter clock sense)
    :param m: the position of the bit of the tile AFTER transformation is performed (k,l)
    :return: the binary value of interest
    """
    # clone list to prevent manipulating the origin data structure
    cur_bit_stream = list(hcq_tiles[i])
    if k:
        cur_bit_stream = apply_swap_operation(cur_bit_stream)
    if l:
        cur_bit_stream = apply_rotation(cur_bit_stream, l)
    return cur_bit_stream[m]


def do_preprocessing(hcq_tiles, I, K, L, M):
    """
    calculate the 4d-param P_iklm to represent the binary patterns of the edges of particles based on their position
    and orientation.
    the param P_iklm is exactly 1, if part i turns out to have a 1-binary value at position m, if it was swapped (k=1)
    or not (k=0) and rotated counter clock wise by l-times 90 degrees compared to its origin position as it was
    provided in the problem instance.

    :param hcq_tiles: a list of lists containing the border patterns of each part
    :return: the 4d-param P_iklm
    """

    ret_val = np.zeros((len(I), len(K), len(L), len(M)), dtype=int)
    for i in I:
        for k in K:
            for l in L:
                for m in M:
                    ret_val[i, k, l, m] = get_binary_value_from_stream(hcq_tiles, i, k, l, m)
    return ret_val


# *******************************************************************************************************************
# * Testcases to validate if code is working as expected                                                            *
# *******************************************************************************************************************

def test_read_hcq_file():
    """
    test whether the read function works as expected

    :return:
    """
    file_path = './data/HC_BQ_01.hcq'
    hcq_tiles = read_hcq_file(file_path)
    assert len(hcq_tiles), "Expected hcq_tiles to contain 6 lists, but found else."
    # make sure, tour distance summary equals to the 4x5 units pieces + 4x5xsqrt(2) units pieces
    assert hcq_tiles[5] == [0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0], \
        f"Expected hcq_tiles[5] to be [0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0], but found else."


def test_apply_swap_operation():
    """
    test, whether the swap operation swaps correctly (horizontal flip of the given part)
    :return:
    """
    cur_bit_stream = [0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0]
    exp_bit_stream = [0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0]
    cur_bit_stream = apply_swap_operation(cur_bit_stream)
    assert cur_bit_stream == exp_bit_stream, "expected bit stream to be " + str(exp_bit_stream) + " but failed."
    # 2nd swap returns origin solution
    exp_bit_stream = [0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0]
    cur_bit_stream = apply_swap_operation(cur_bit_stream)
    assert cur_bit_stream == exp_bit_stream, "expected bit stream to be " + str(exp_bit_stream) + " but failed."

    cur_bit_stream = [1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0]
    exp_bit_stream = [0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0]
    cur_bit_stream = apply_swap_operation(cur_bit_stream)
    assert cur_bit_stream == exp_bit_stream, "expected bit stream to be " + str(exp_bit_stream) + " but failed."
    # 2nd swap returns origin solution
    exp_bit_stream = [1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0]
    cur_bit_stream = apply_swap_operation(cur_bit_stream)
    assert cur_bit_stream == exp_bit_stream, "expected bit stream to be " + str(exp_bit_stream) + " but failed."


def test_apply_rotation():
    """
    test whether the rotation operation works as expected and turns the tiles in l steps of 90 degrees counter
     clock wise manner
    :return:
    """
    cur_bit_stream = [1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0]
    exp_bit_stream = [0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0]
    cur_bit_stream = apply_rotation(cur_bit_stream, 1)
    assert cur_bit_stream == exp_bit_stream, "expected bit stream to be " + str(exp_bit_stream) + " but failed."
    cur_bit_stream = [1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0]
    exp_bit_stream = [1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1]
    cur_bit_stream = apply_rotation(cur_bit_stream, 2)
    assert cur_bit_stream == exp_bit_stream, "expected bit stream to be " + str(exp_bit_stream) + " but failed."
    cur_bit_stream = [1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0]
    exp_bit_stream = [0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0]
    cur_bit_stream = apply_rotation(cur_bit_stream, 3)
    assert cur_bit_stream == exp_bit_stream, "expected bit stream to be " + str(exp_bit_stream) + " but failed."
    # turn tile 4 times 90 degrees returns origin stream
    cur_bit_stream = [1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0]
    exp_bit_stream = [1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0]
    cur_bit_stream = apply_rotation(cur_bit_stream, 4)
    assert cur_bit_stream == exp_bit_stream, "expected bit stream to be " + str(exp_bit_stream) + " but failed."


def test_do_preprocessing():
    file_path = './data/HC_BQ_01.hcq'
    hcq_tiles = read_hcq_file(file_path)
    # define index sets
    I = range(len(hcq_tiles))  # the ID (1-6) of the tiles as they are currently placed on the table
    # 1, 2, 5, 6 (vertically), 3, 2, 4 (horizontally)
    K = range(2)  # apply horizontal swap to tile before placing in the net (=1) or not (=0)?
    L = range(5)  # rotate tile counterclockwise for 0, 1, 2, 3, 4 x 90° ?
    M = range(20)  # bit-string length for every part (4*5)

    # define parameters
    P = do_preprocessing(hcq_tiles, I, K, L, M)

    # take part 1, do not swap and do not turn, expect origin line 1 in file
    exp_bit_stream = [0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0]
    assert P[0, 0, 0, :].tolist() == exp_bit_stream
    # take part 1, swap and do not turn, expect swapped bit-string
    exp_bit_stream = [0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0]
    assert P[0, 1, 0, :].tolist() == exp_bit_stream
    # take part 1, do not swap and turn 1x cclk wise, expect rotated bit-string
    exp_bit_stream = [0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0]
    assert P[0, 0, 1, :].tolist() == exp_bit_stream
    # take part 5, do not swap and do not turn cclk wise, expect rotated bit-string
    exp_bit_stream = [0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0]
    assert P[4, 0, 0, :].tolist() == exp_bit_stream


# *******************************************************************************************************************
# * Main Program                                                                                                    *
# *******************************************************************************************************************

def main():
    file_paths = ['./data/HC_BQ_01.hcq', './data/HC_BQ_02.hcq', './data/HC_BQ_03.hcq']

    for file in file_paths:
        instance_name = os.path.splitext(os.path.basename(file))[0]
        hcq_tiles = read_hcq_file(file)
        print(f"*** Solving File {file} with MILP:")
        solution = solve_hcq_milp(hcq_tiles)
        instructions = [solution[key] for key in sorted(solution)]
        print("Solution:")
        for instruction in instructions:
            print(f"- {instruction}")


if __name__ == "__main__":
    main()

The whole script is also stored here on my github account.

Validate the Model

The best way to validate the model and implementation is by solving a few test instances. The table below shows some of these test cases. I’ve successfully applied the resulting instructions to the different instances multiple times, and the cubes were successfully solved each time.

Test Instance 1: HC_BQ_001 (see file here)

Test Instance Description

Initial setup of test instance **HC_BQ_001**

Soltuion Instructions:
Take part 6 and put it on place 1.
Take part 5 rotate it 1 times 90 degrees counter-clockwise and put it on place 2.
Take part 3 rotate it 3 times 90 degrees counter-clockwise and put it on place 3.
Take part 2 rotate it 3 times 90 degrees counter-clockwise and put it on place 4.
Take part 1 rotate it 3 times 90 degrees counter-clockwise and put it on place 5.
Take part 4 and put it on place 6.

*How to apply the solver instructions to assemble the cube from instance HC_BQ_001*

Once the pieces are correctly placed, the surfaces can be folded upwards so that the piece in position 2 forms the base, and the piece from position 6 becomes the top surface.

Test Instance 2: HC_BQ_002 (see file here)

Soltuion Instructions:
Take part 6 and put it on place 1.
Take part 4 rotate it 3 times 90 degrees counter-clockwise and put it on place 2.
Take part 5 rotate it 3 times 90 degrees counter-clockwise and put it on place 3.
Take part 1 and put it on place 4.
Take part 2 rotate it 1 times 90 degrees counter-clockwise and put it on place 5.
Take part 3 and put it on place 6

Test Instance 3: HC_BQ_003 (see file here)
Same as instance 2, but other scrambling of tiles:

Solution Instructions:
Take part 6 rotate it 1 times 90 degrees counter-clockwise and put it on place 1.
Take part 4 rotate it 1 times 90 degrees counter-clockwise and put it on place 2.
Take part 2, swap it horizontally, rotate it 1 times 90 degrees counter-clockwise and put it on place 3.
Take part 5, swap it horizontally, rotate it 3 times 90 degrees counter-clockwise and put it on place 4.
Take part 3 rotate it 3 times 90 degrees counter-clockwise and put it on place 5.
Take part 1 and put it on place 6.

Conclusion

Whether it was worth writing the solver as a mathematical model is debatable, especially since the source code with the mathematical model is larger than the one for the micro:bit, where the problem was solved using dynamic programming and backtracking.

That said, it was fun, and the idea of a mathematical model solving 3D puzzles is pretty cool. 🙂