-
-
Save phabee/acc312711c5b2d3108cfe2e5d9faa19b to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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