Skip to content

Instantly share code, notes, and snippets.

@phabee
Last active September 1, 2024 18:27
Show Gist options
  • Save phabee/acc312711c5b2d3108cfe2e5d9faa19b to your computer and use it in GitHub Desktop.
Save phabee/acc312711c5b2d3108cfe2e5d9faa19b to your computer and use it in GitHub Desktop.
# 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment