```
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import ipywidgets as widgets
from IPython.display import HTML
from math import ceil
```

The Soma cube is the brainchild of Danish mathematician, inventor, designer, writer and poet Piet Hein. The puzzle consists of seven pieces and the goal is to assemble these pieces into a 3x3x3 cube or other shapes.

Piet Hein came up with this puzzle during a lecture on quantum physics by Werner Heisenberg. Instead of paying attention to the renowned theoretical physicist, Piet prefered idling away his time by thinking up frivolous brainteasers. Such a waste, not to mention rude! What’s even worse is that people seemed to actually enjoy his puzzle. This meant that Piet Hein’s invention wasn’t just a time-waster for him but for everyone who got hooked on it.

Initially, the Soma cube was mostly known in Scandinavian countries, but things took a dark turn when Martin Gardner featured it in his Mathematical Games column in Scientific American in September 1958. Suddenly, the whole world got introduced to this time-waster.

As a self-proclaimed restorer of productivity, I hate to see people trifle away their lives trying to solve these inane puzzles. For this reason I’ve created a Soma solver in Python so people can go back to spending time on more fruitful pursuits.

Read on to learn how I did this.

## The Soma puzzle pieces

Each of the seven pieces has a different shape. I’ve given each shape a different color to ease visualizing solutions later on.

As mentioned above, the objective is to assemble these pieces into various shapes, like a cube.

## Importing libraries

Let’s start by importing some libraries. We’ll use `matplotlib`

to visualize solutions in 3D.

## Representing the pieces

We’ll represent each piece as a list of coordinates in 3D space. Each tuple is an (x, y, z) coordinate.

```
= [(0,0,0),(1,0,0),(1,1,0),(2,1,0)] # blue piece
z = [(0,0,0),(0,1,0),(0,1,1),(1,1,0)] # red piece
p = [(0,0,0),(0,1,0),(1,1,0),(0,2,0)] # purple piece
t = [(0,0,0),(1,0,0),(0,1,0),(0,1,1)] # brown piece
b = [(0,0,0),(0,0,1),(0,1,0),(1,1,0)] # yellow piece
a = [(0,0,0),(1,0,0),(2,0,0),(0,1,0)] # orange piece
l = [(0,0,0),(1,0,0),(0,1,0)] # green piece v
```

The letters `z`

, `t`

, `p`

, `b`

, `a`

, `l`

, and `v`

look (with some imagination) like the pieces.

Let’s put all pieces and colors in a list so we can access them by index.

```
= [z, p, t, b, a, l, v]
pieces = ["blue", "red", "purple", "brown", "yellow", "orange", "green"] colors
```

## Visualizing the pieces and solution

Before writing the solver, let’s consider how we can create a 3D visualization of the Soma cube. I’ve chosen to use Matplotlib for its 3D rendering capabilities.

Matplotlib can draw voxels (3D pixels) which we can use to visualize a single piece, multiple pieces, or even the full solution in 3D. The following function takes a 3D numpy array that represents the voxels and plots them.

```
def plot_solution(colors, ax=None):
if not ax:
= plt.figure()
fig
# axis with 3D projection
= fig.add_subplot(projection='3d')
ax
'equal')
ax.set_aspect(
ax.set_axis_off()
# draw each voxel with a color (each voxel unequal to None)
= (colors != None)
voxels
=colors, edgecolors=colors) ax.voxels(voxels, facecolors
```

To use this function we need to call it with a 3D numpy array that represents each voxel.

`= np.empty((3,3,3), dtype='object') voxels `

Each element in this array is initialized to `None`

.

` voxels`

```
array([[[None, None, None],
[None, None, None],
[None, None, None]],
[[None, None, None],
[None, None, None],
[None, None, None]],
[[None, None, None],
[None, None, None],
[None, None, None]]], dtype=object)
```

Let’s set each element of `voxels`

to a color and call `plot_solution`

. We can set elements seperately, like this:

```
0][0][0] = 'yellow'
voxels[1][0][0] = 'yellow'
voxels[0][0][1] = 'yellow'
voxels[0][1][1] = 'yellow' voxels[
```

Or we can set all elements at once.

```
= np.array([
voxels 'yellow', 'yellow', 'orange'],
[['brown', 'yellow', 'orange'],
['brown', 'brown', 'orange']],
['yellow', 'green', 'orange'],
[['brown', 'blue', 'blue'],
['blue', 'blue', 'red']],
['purple', 'green', 'green'],
[['purple', 'purple', 'red'],
['purple', 'red', 'red']]
[ ])
```

```
plot_solution(voxels) plt.show()
```

This is the same solution as shown above and one of the 240 possible unique ways of packing the seven Soma pieces into a 3x3x3 cube.

Let’s see how we can visualize a single piece.

```
def plot_piece(piece, color, ax=None):
= np.max([np.max(piece) + 1, 3])
max_dim = np.empty((max_dim, max_dim, max_dim), dtype='object')
voxels for x, y, z in piece:
= color
voxels[x][y][z] plot_solution(voxels, ax)
```

```
"blue")
plot_piece(z, plt.show()
```

## Animating a solution

In the solution above it was impossible to tell how some pieces were positioned. The following function animates the assembly of a solution by adding pieces one by one.

```
def animate_solution(colors):
= colors.shape
w, d, h
= []
show_order for a in range(1,max(w, d, h)+1):
for z in range(min(a,h)):
for x in range(min(a,w)):
for y in range(d-1,d-1-min(a,d),-1):
= colors[x][y][z]
color if not (color is None or color in show_order):
show_order.append(color)
= plt.figure()
fig
= fig.add_subplot(projection='3d')
ax
def update(frame):
ax.clear()'equal')
ax.set_aspect(
ax.set_axis_off()
= np.in1d(colors, show_order[:frame+1]).reshape(colors.shape)
voxels
=colors, edgecolors=colors)
ax.voxels(voxels, facecolors
plt.close()
= animation.FuncAnimation(fig, update, frames=len(show_order), interval=1000)
anim
return anim
```

```
= animate_solution(voxels)
anim ='once')) HTML(anim.to_jshtml(default_mode
```

Now it’s much clearer how a solution is constructed.

## Rotating pieces

We will implement the solver using a simple recursive backtracking algorithm. The solver will try to fit the pieces in their various orientations. For this reason, we define some helper functions to rotate a piece around an axis. I also defined an `identity`

function that just returns the input as is. Its use will become clear later.

```
= lambda cubelets: [(x, z, -y) for (x, y, z) in cubelets]
rotate_x = lambda cubelets: [(z, y, -x) for (x, y, z) in cubelets]
rotate_y = lambda cubelets: [(-y, x, z) for (x, y, z) in cubelets]
rotate_z = lambda cubelets: cubelets identity
```

For instance, to rotate piece `z`

around the x-axis we can make use of `rotate_x`

.

` rotate_x(z)`

`[(0, 0, 0), (1, 0, 0), (1, 0, -1), (2, 0, -1)]`

We can see that rotating a piece may make some coordinates negative. The following function translates a piece so that all `x`

, `y`

, `z`

coordinates are minimal but positive.

```
def translate(piece):
= np.min(np.array(piece), axis=0) * -1
d_x, d_y, d_z return [(x + d_x, y + d_y, z + d_z) for (x, y, z) in piece]
```

Now we can do the same rotation as above but end up with coordinates that are all positive.

` translate(rotate_x(z))`

`[(0, 0, 1), (1, 0, 1), (1, 0, 0), (2, 0, 0)]`

We will use this `translate`

function when we display individual pieces.

## Generating all orientations for each piece

Using the functions defined above we can generate all orientations for each piece. Several transformations are performed one after another. The `identity`

function just returns the original orientation.

```
def generate_rotations(piece):
= []
orientations for f_a in [identity, rotate_x, rotate_y, rotate_z]:
for f_b in [identity, rotate_x, rotate_y, rotate_z]:
for f_c in [identity, rotate_x, rotate_y, rotate_z]:
for f_d in [identity, rotate_x, rotate_y, rotate_z]:
for f_e in [identity, rotate_x, rotate_y, rotate_z]:
= sorted(f_a(f_b(f_c(f_d(f_e(piece))))))
rot_piece = rot_piece[0]
min_x, min_y, min_z = [(x - min_x, y - min_y, z - min_z) for x, y, z in rot_piece]
trans_rot_piece if trans_rot_piece not in orientations:
orientations.append(trans_rot_piece)
return orientations
```

We can call this function for one piece.

` generate_rotations(z)`

```
[[(0, 0, 0), (1, 0, 0), (1, 1, 0), (2, 1, 0)],
[(0, 0, 0), (1, 0, -1), (1, 0, 0), (2, 0, -1)],
[(0, 0, 0), (0, 0, 1), (0, 1, -1), (0, 1, 0)],
[(0, 0, 0), (0, 1, 0), (1, -1, 0), (1, 0, 0)],
[(0, 0, 0), (1, -1, 0), (1, 0, 0), (2, -1, 0)],
[(0, 0, 0), (0, 1, 0), (0, 1, 1), (0, 2, 1)],
[(0, 0, 0), (0, 0, 1), (1, 0, 1), (1, 0, 2)],
[(0, 0, 0), (0, 1, -1), (0, 1, 0), (0, 2, -1)],
[(0, 0, 0), (1, 0, 0), (1, 0, 1), (2, 0, 1)],
[(0, 0, 0), (0, 1, 0), (1, 1, 0), (1, 2, 0)],
[(0, 0, 0), (0, 0, 1), (0, 1, 1), (0, 1, 2)],
[(0, 0, 0), (0, 0, 1), (1, 0, -1), (1, 0, 0)]]
```

Or we can apply this function to every piece in the `pieces`

list we defined above.

`= list(map(generate_rotations, pieces)) orientations `

We can display all orientations for each piece. Note how we apply `translate`

to make sure all coordinates end up in the positive quadrant.

```
= 8
cols = sum(list(map(len,orientations)))
total_number_of_orientations = ceil(total_number_of_orientations / cols)
rows = plt.subplots(rows, cols, figsize=(cols*3,rows*3),
fig, axs = { 'projection':'3d' }, gridspec_kw = {'wspace':0.0, 'hspace':0.0})
subplot_kw
=0.96)
fig.subplots_adjust(topf"Orientations", fontsize=16)
fig.suptitle(
= 0
count for i in range(len(pieces)):
for j in range(len(orientations[i])):
// cols, count % cols])
plot_piece(translate(orientations[i][j]), colors[i], axs[count += 1
count
while count < cols * rows:
// cols, count % cols].set_axis_off()
axs[count += 1
count
plt.show()
```

In the solver below, for each piece we’ll see which orientation fits where.

## Writing the solver

We’re finally ready to work on the solver for the Soma cube. The solver will take in a list of coordinates (that make up the cube that we want to fill up with pieces) and recursively search for a solution.

We can generate this list of coordinates using a simple list comprehension.

`= [(x, y, z) for x in range(3) for y in range(3) for z in range(3)] cube_coordinates `

` cube_coordinates`

```
[(0, 0, 0),
(0, 0, 1),
(0, 0, 2),
(0, 1, 0),
(0, 1, 1),
(0, 1, 2),
(0, 2, 0),
(0, 2, 1),
(0, 2, 2),
(1, 0, 0),
(1, 0, 1),
(1, 0, 2),
(1, 1, 0),
(1, 1, 1),
(1, 1, 2),
(1, 2, 0),
(1, 2, 1),
(1, 2, 2),
(2, 0, 0),
(2, 0, 1),
(2, 0, 2),
(2, 1, 0),
(2, 1, 1),
(2, 1, 2),
(2, 2, 0),
(2, 2, 1),
(2, 2, 2)]
```

The solver will yield the same list but each tuple will have an extra element indicating the color that coordinate will get after all pieces have been placed.

We’ll use a list of booleans to keep track of which pieces have been placed.

`= [False] * 7 pieces_used `

```
def solve_soma_(solution, i):
if i == 27:
yield solution
else:
= solution[i]
x, y, z, _ for piece in range(7):
if not pieces_used[piece]:
for orientation in orientations[piece]:
= [(x + d_x, y + d_y, z + d_z, None) for (d_x, d_y, d_z) in orientation]
empty_coords if all([tup in solution for tup in empty_coords]):
# placing piece: replace None with color
= True
pieces_used[piece] = [(x + d_x, y + d_y, z + d_z, colors[piece]) for (d_x, d_y, d_z) in orientation]
filled_coords = sorted([tup for tup in solution if tup not in empty_coords] + filled_coords)
new_solution
# find next empty coordinate
= i
j while j < 27 and new_solution[j][3]:
+= 1
j
# continue search
yield from solve_soma_(new_solution, j)
= False pieces_used[piece]
```

```
def solve_soma(coordinates):
global pieces_used
= [False] * 7
pieces_used
= sorted([(x, y, z, None) for x, y, z in coordinates])
solution yield from solve_soma_(solution, 0)
```

You call the `solve_soma`

function with a list of coordinates representing the space that will be packed with Soma pieces. The actual algorithm is implemented in `solve_soma_`

. The algorithm does the following:

- It checks if all coordinates are filled. If so, it yields the solution.
- For each possible orientation of a piece, it checks if it can be placed from the current coordinate without overlapping with filled coordinates that are already in the solution.
- If a valid placement is found, it recursively explores how to fill from the next empty coordinate.
- After exploring a piece and its orientation, it tries other possibilities.

Different solutions may be duplicates of each other in the sense that one solution is a rotation of another.

## Finding and displaying a solution

The solver is a generator function that yields all solutions one by one. Let’s start by just displaying the first solution it finds.

The following function finds the first solution given a list of coordinates to fill and plots this solution.

```
def solve_and_plot_first_solution(coordinates, ax=None):
# instantiate generator
= solve_soma(coordinates)
gen
# find first solution
= next(gen)
solution
# find voxel space size
= np.max(coordinates) + 1
max_dim
# fill 3D voxel array and plot it
= np.empty((max_dim, max_dim, max_dim), dtype='object')
voxels for x, y, z, color in solution:
= color
voxels[x][y][z]
plot_solution(voxels, ax)
```

```
solve_and_plot_first_solution(cube_coordinates) plt.show()
```

We’ll also define a function that finds the first solution and then animates it.

```
def solve_and_animate_first_solution(coordinates):
# instantiate generator
= solve_soma(coordinates)
gen
# find first solution
= next(gen)
solution
# find voxel space size
= np.max(coordinates) + 1
max_dim
# fill 3D voxel array and plot it
= np.empty((max_dim, max_dim, max_dim), dtype='object')
voxels for x, y, z, color in solution:
= color
voxels[x][y][z]
= animate_solution(voxels)
anim ='once'))) display(HTML(anim.to_jshtml(default_mode
```

` solve_and_animate_first_solution(cube_coordinates)`

## Creating other figures

We called the solver with `cube_coordinates`

that contains the coordinates that make up a cube. Instead of a cube, we can also pass a sorted list of coordinates that make up another shape.

Let’s look at some examples.

```
= sorted(
pyramid_coordinates 0) for x in range(1, 4) for y in range(5)] +
[(x, y, 0) for x in [0, 4] for y in range(1, 4)] +
[(x, y, 2, y, 1) for y in range(1, 4)] +
[(1, 2, 1), (3, 2, 1), (2, 2, 2)]
[(
)
= sorted(
turtle_coordinates for x in range(1, 4) for y in range(1, 4) for z in range(2)] +
[(x, y, z) 4, 2, z) for z in range(3)] +
[(0, 2, 0), (1, 0, 0), (1, 4, 0), (3, 0, 0), (3, 4, 0), (5, 2, 2)]
[(
)
= sorted(
tower_coordinates for x in range(2) for y in range(2) for z in range(7) if (x, y, z) != (1, 0, 6)]
[(x, y, z)
)
= sorted(
bear_coordinates 0) for x in range(3) for y in range(2)] +
[(x, y, 1, z) for x in range(3) for z in range(1, 6)] +
[(x, 0, z) for x in [0, 2] for z in [1, 3, 4]]
[(x,
)
= sorted(
tunnel_coordinates 0) for x in range(5) for y in range(3) if x != 2] +
[(x, y, 1) for x in range(1, 4) for y in range(3) if x != 2] +
[(x, y, 2) for x in range(1, 4) for y in range(3)]
[(x, y,
)
= sorted(
tub_coordinates 0) for x in range(5) for y in range(3)] +
[(x, y, 1) for x in range(5) for y in range(3) if x in [0, 4] or y in [0, 2]]
[(x, y,
)
= sorted(
airplane_coordinates 0) for x in range(1, 6) for y in range(7) if x == 3 or y in [0, 1]] +
[(x, y, 1) for x in range(7) for y in range(7) if y == 0 or (x == 3 and y != 5)]
[(x, y,
)
= sorted(
bench_coordinates for x in range(5) for y in range(2) for z in range(3) if y == 1 or x in [0, 4]] +
[(x, y, z) 0, 1) for x in range(1, 4)] +
[(x, 1, 3) for x in range(1, 4)]
[(x,
)
= sorted(
duck_coordinates for x in range(4) for y in range(3) for z in range(2)] +
[(x, y, z) 3, 1, 2), (3, 1, 3), (4, 1, 3)]
[(
)
= sorted(
cascade_coordinates for x in range(3) for y in range(3) for z in range(5) if (2-x) + y >= z]
[(x, y, z)
)
= sorted(
chair_coordinates for x in range(3) for y in range(3) for z in range(2)] +
[(x, y, z) 2, z) for x in range(3) for z in range(2, 5)]
[(x,
)
= sorted(
castle_coordinates 0) for x in range(5) for y in range(5) if (x, y) not in [(0,2), (4,2)]] +
[(x, y, 0, 0, 1), (0, 4, 1), (4, 0, 1), (4, 4, 1)]
[(
)
= sorted(
dog_coordinates 0) for x in range(6) for y in range(3) if (x, y) not in [(3, 0), (3, 2), (5, 1)]] +
[(x, y, 1) for x in range(5) for y in range(3) if x == 4 or y == 1] +
[(x, y, 1, z) for x in range(3, 6) for z in range(2, 4) if (x, z) != (5, 3)]
[(x,
)
= sorted(
cross_coordinates for x in range(3) for y in range(3) for z in range(2)] +
[(x, y, z) 1, y, 2) for y in range(3)] +
[(1, 1, z) for z in range(3, 7)] +
[(0, 1, 5), (2, 1, 5)]
[(
)
= sorted(
elephant_coordinates 0, 0, 0), (3, 0, 0), (0, 2, 0), (3, 2, 0)] +
[(1) for x in range(4) for y in range(3) if (x, y) != (3, 1)] +
[(x, y, 1, z) for x in range(5) for z in range(2, 4) if (x, z) != (0, 3)] +
[(x, 4, 1, 1), (2, 0, 2), (2, 2, 2)]
[(
)
= sorted(
stairs_coordinates for y in range(3) for z in range(3) for x in range(4-z)]
[(x, y, z)
)
= sorted(
snake_coordinates for y in range(4) for z in range(2) for x in range(2*(3-y), (2*(3-y))+3)] +
[(x, y, z) 0, z) for x in range(8, 10) for z in range(2, 4) if x == 8 or z == 3]
[(x,
)
= sorted(
skyscraper_coordinates 0) for x in range(3) for y in range(3)] +
[(x, y, for x in range(2) for y in range(1, 3) for z in range(1, 5)] +
[(x, y, z) 1, 1, 5), (1, 1, 6)]
[(
)
= sorted(
wall_coordinates - 1, 0) for x in range(1, 6)] +
[(x, x for x in range(6) for z in range(2)] +
[(x, x, z) + 1, z) for x in range(5) for z in range(2)]
[(x, x )
```

Let’s place these figures and their names into lists.

```
= ["Cube", "Pyramid", "Turtle", "Tower", "Bear", "Tunnel", "Tub", "Airplane",
figure_names "Bench", "Duck", "Cascade", "Chair", "Castle", "Dog", "Cross", "Elephant",
"Stairs", "Snake", "Skyscraper", "Wall"
]= [
figure_coordinates
cube_coordinates, pyramid_coordinates, turtle_coordinates,
tower_coordinates, bear_coordinates, tunnel_coordinates,
tub_coordinates, airplane_coordinates, bench_coordinates,
duck_coordinates, cascade_coordinates, chair_coordinates,
castle_coordinates, dog_coordinates, cross_coordinates,
elephant_coordinates, stairs_coordinates, snake_coordinates,
skyscraper_coordinates, wall_coordinates ]
```

You can add more figures if you want. Just make sure that each list contains exactly 27 coordinates.

We can plot all figures as subplots.

```
= 5
cols = len(figure_coordinates)
number_of_figures = ceil(number_of_figures / cols)
rows = plt.subplots(rows, cols, figsize=(cols*3,rows*3),
fig, axs = { 'projection':'3d' }, gridspec_kw = {'wspace':0.5, 'hspace':0.0})
subplot_kw
f"Figures", fontsize=16)
fig.suptitle(
= 0
count for i in range(number_of_figures):
// cols, count % cols])
solve_and_plot_first_solution(figure_coordinates[i], axs[count // cols, count % cols].set_title(figure_names[i])
axs[count += 1
count
while count < cols * rows:
// cols, count % cols].set_axis_off()
axs[count += 1
count
plt.show()
```

I have to admit that some figures are pretty cute.

## Showing solutions one by one

Most figures can be assembled in more than one way. Using some widgets we can put together a simple GUI that lets us select a figure, solve it and iterate over the solutions.

We’ll use a dropdown menu, a button, an output widget, and a label widget. The solutions are shown as an animation we can step through.

```
= widgets.Dropdown(options=zip(figure_names, figure_coordinates), value=cube_coordinates, description="Figure: ")
figure_dropdown = widgets.Button(description="Next")
next_button = widgets.HBox([figure_dropdown, next_button])
hbox = widgets.Output()
output = widgets.Label()
label
display(hbox, output, label)
```

To keep track of which solution to show next, we declare the following global variables:

```
= None
gen = None solution
```

The gen is a generator object that generates solutions for a given figure using the `solve_soma`

function. The solution variable stores the next solution to be displayed. It is updated each time the “Next” button is clicked to display the next solution.

```
def figure_dropdown_change(change):
global gen
global solution
if change['type'] == 'change' and change['name'] == 'value':
= change['new']
value = 'none'
next_button.layout.display
output.clear_output()
= solve_soma(value)
gen try:
= next(gen)
solution with output:
# find voxel space size
= np.max(value) + 1
max_dim
# fill 3D voxel array and plot it
= np.empty((max_dim, max_dim, max_dim), dtype='object')
voxels for x, y, z, color in solution:
= color
voxels[x][y][z]
= animate_solution(voxels)
anim ='once')))
display(HTML(anim.to_jshtml(default_modeexcept StopIteration:
= "No solutions found."
label.value return
try:
= next(gen).copy()
solution = ""
label.value = 'inline-block'
next_button.layout.display except StopIteration:
= ""
label.value
figure_dropdown.observe(figure_dropdown_change)
```

The `figure_dropdown_change`

function handles changes to the figure dropdown menu. It clears the output widget, displays the first solution, and stores the next solution in the solution variable (if any).

```
def on_next_button_clicked(b):
global gen
global solution
=True)
output.clear_output(waitwith output:
# find voxel space size
= np.max(figure_dropdown.value) + 1
max_dim
# fill 3D voxel array and plot it
= np.empty((max_dim, max_dim, max_dim), dtype='object')
voxels for x, y, z, color in solution:
= color
voxels[x][y][z]
= animate_solution(voxels)
anim ='once')))
display(HTML(anim.to_jshtml(default_mode
try:
= next(gen)
solution except StopIteration:
= ""
label.value = 'none'
next_button.layout.display
next_button.on_click(on_next_button_clicked)
'type': 'change', 'name': 'value', 'new': figure_dropdown.value}) figure_dropdown_change({
```

The `on_next_button_clicked`

function handles clicks on the “Next” button by clearing the output widget, displaying the current solution, and generating the next solution.

## Solving your own Soma puzzles

There you have it! No more wasting time on solving Soma puzzles. If you want to experiment with an interactive version of the solver, you can check out on of the following links:

Now go do something useful!