How Many Boxes are on the Trailer?

I recently came across an interesting puzzle:

Without shading or a 3D-view, we can not give a definitive answer. Taking a look at the side view, we can give an upper bound. If we expect, from the side view, every visible cube to have two invisible boxes behind it, we get to 51 cubes. But what is the lower bound?

I thought, it would be a good idea and exercise, to solve this problem using a constraint solver. First I needed to decide on a model. As there exists a grid of $7*3*3$ fields where each contains either a box or it does not. That sounds suspiciously like something Boolean. So let's model it as a Boolean satisfiability problem (SAT).

Since I haven't before, I wanted to try out PySAT and take it for a spin: The SAT solvers work with Boolean variables without names, just numbers. Additionally, our interface will create additional variables when we create more complex formulas due to Tseytin transformation. Tseytin transformation essentially creates a variable representing a subformula and the constraints for equivalence between the valuation of the newly created variable and the subformula. If you now want to create a new formula that contains the original formula as a subformula, you can refer to the single newly created variable, representing the whole subformula.

Because, as I said, the interface with the solver will be numbers-based, we need to decide on some representation or indexing between the solver's variables and our domain. We'll use a namedtuple for a Position P on the trailer, where x is along the depth of the trailer when viewed from the back, y is along the width of the trailer when viewed from the back, and z is along the height of the trailer, with the origin being at the bottom front left. We also need an empty list of constraints that we add to along the way:

from pysat.formula import *
from collections import namedtuple

P = namedtuple('P', ['x', 'y', 'z'])

boxes = {}

depth = 7
width = 3
height = 3

for x in range(depth):
  for y in range(width):
    for z in range(height):
      boxes[(x,y,z)] = Atom(P(x,y,z))

constraints = []

Instead of a namedtuple i initially wanted to use anonymous 3-tuples, but a bug in the library prevented me from doing that.

First we want to create the constraints for the side view. If we can see a box from that side, we know that at least one space in that group of 3 is filled by a box. We can simply use a disjunction over those three variables. If the space is empty, we know that there are no boxes. We can negate the disjunction. We "draw" the shape to make sure we got it right:

# side view
X = True
O = False
side_view = [
  [X,X,X,X,O,O,O],
  [X,X,X,X,X,X,O],
  [X,X,X,X,X,X,X],
]
for x in range(depth):
  for z in range(height):
    if side_view[-z-1][x]:
      constraints.append(Or(*[boxes[(x,y,z)] for y in range(width)]))
    else:
      constraints.append(Neg(Or(*[boxes[(x,y,z)] for y in range(width)])))

For the top view and back view we do not need to draw the shape that we are seeing, because every row/column/stack has a box in it and we do not need to consider the other case.

# back view
for y in range(width):
  for z in range(height):
    constraints.append(Or(*[boxes[(x,y,z)] for x in range(depth)]))

# top view
for y in range(width):
  for x in range(depth):
    constraints.append(Or(*[boxes[(x,y,z)] for z in range(height)]))

This looks like we are finished, but we have not yet considered gravity. If a box is anywhere not on the bottom level, it has to have another box under it. That is an implication:

# gravity
for x in range(depth):
  for y in range(width):
    for z in range(height-1):
      constraints.append(Implies(boxes[(x,y,z+1)], boxes[(x,y,z)]))

When I tried this out initially, there was a bug in the pySAT library where you could not use implications that are not on the top-level. This has been fixed quite fast.

Now we want to know the minimum number of needed boxes. We will get one solution, check how many boxes we needed, and add a constraint that we want one box less. We will do this until our formula is no longer satisfiable.

f = And(*constraints)

from pysat.solvers import Solver
with Solver(bootstrap_with=f) as solver:
  while solver.solve():
    model = solver.get_model()
    atoms = Formula.formulas((b for b in model if b > 0), atoms_only=True)
    boxes_needed = len(atoms)
    print("Solution with", boxes_needed, "boxes found")
    print(atoms)
    print("="*80)

    from pysat.card import CardEnc
    solver.append_formula(CardEnc.atmost(f.literals(boxes.values()), boxes_needed-1, vpool=Formula.export_vpool()))

It is important, when counting boxes, that we only extract the formulas/atoms that we created and not the ones created internally due to Tseyitin transformation. We count them and add a cardinality constraint, where we say, that we want at most one box less than we currently have used. An after only three iterations we get the following result:

Solution with 33 boxes found
[Atom(P(x=0, y=0, z=0)), Atom(P(x=0, y=1, z=0)), Atom(P(x=0, y=2, z=0)), Atom(P(x=0, y=1, z=1)), Atom(P(x=0, y=2, z=1)), Atom(P(x=0, y=1, z=2)), Atom(P(x=0, y=2, z=2)), Atom(P(x=1, y=0, z=0)), Atom(P(x=1, y=1, z=0)), Atom(P(x=1, y=2, z=0)), Atom(P(x=1, y=0, z=1)), Atom(P(x=1, y=0, z=2)), Atom(P(x=2, y=0, z=0)), Atom(P(x=2, y=1, z=0)), Atom(P(x=2, y=2, z=0)), Atom(P(x=2, y=0, z=1)), Atom(P(x=2, y=0, z=2)), Atom(P(x=3, y=0, z=0)), Atom(P(x=3, y=1, z=0)), Atom(P(x=3, y=2, z=0)), Atom(P(x=3, y=0, z=1)), Atom(P(x=3, y=0, z=2)), Atom(P(x=4, y=0, z=0)), Atom(P(x=4, y=1, z=0)), Atom(P(x=4, y=2, z=0)), Atom(P(x=4, y=0, z=1)), Atom(P(x=5, y=0, z=0)), Atom(P(x=5, y=1, z=0)), Atom(P(x=5, y=2, z=0)), Atom(P(x=5, y=0, z=1)), Atom(P(x=6, y=0, z=0)), Atom(P(x=6, y=1, z=0)), Atom(P(x=6, y=2, z=0))]
================================================================================
Solution with 32 boxes found
[Atom(P(x=0, y=0, z=0)), Atom(P(x=0, y=1, z=0)), Atom(P(x=0, y=2, z=0)), Atom(P(x=0, y=0, z=1)), Atom(P(x=0, y=0, z=2)), Atom(P(x=1, y=0, z=0)), Atom(P(x=1, y=1, z=0)), Atom(P(x=1, y=2, z=0)), Atom(P(x=1, y=1, z=1)), Atom(P(x=1, y=1, z=2)), Atom(P(x=2, y=0, z=0)), Atom(P(x=2, y=1, z=0)), Atom(P(x=2, y=2, z=0)), Atom(P(x=2, y=0, z=1)), Atom(P(x=2, y=2, z=1)), Atom(P(x=2, y=2, z=2)), Atom(P(x=3, y=0, z=0)), Atom(P(x=3, y=1, z=0)), Atom(P(x=3, y=2, z=0)), Atom(P(x=3, y=0, z=1)), Atom(P(x=3, y=0, z=2)), Atom(P(x=4, y=0, z=0)), Atom(P(x=4, y=1, z=0)), Atom(P(x=4, y=2, z=0)), Atom(P(x=4, y=2, z=1)), Atom(P(x=5, y=0, z=0)), Atom(P(x=5, y=1, z=0)), Atom(P(x=5, y=2, z=0)), Atom(P(x=5, y=2, z=1)), Atom(P(x=6, y=0, z=0)), Atom(P(x=6, y=1, z=0)), Atom(P(x=6, y=2, z=0))]
================================================================================
Solution with 31 boxes found
[Atom(P(x=0, y=0, z=0)), Atom(P(x=0, y=1, z=0)), Atom(P(x=0, y=2, z=0)), Atom(P(x=0, y=0, z=1)), Atom(P(x=0, y=0, z=2)), Atom(P(x=1, y=0, z=0)), Atom(P(x=1, y=1, z=0)), Atom(P(x=1, y=2, z=0)), Atom(P(x=1, y=2, z=1)), Atom(P(x=1, y=2, z=2)), Atom(P(x=2, y=0, z=0)), Atom(P(x=2, y=1, z=0)), Atom(P(x=2, y=2, z=0)), Atom(P(x=2, y=0, z=1)), Atom(P(x=2, y=0, z=2)), Atom(P(x=3, y=0, z=0)), Atom(P(x=3, y=1, z=0)), Atom(P(x=3, y=2, z=0)), Atom(P(x=3, y=1, z=1)), Atom(P(x=3, y=1, z=2)), Atom(P(x=4, y=0, z=0)), Atom(P(x=4, y=1, z=0)), Atom(P(x=4, y=2, z=0)), Atom(P(x=4, y=2, z=1)), Atom(P(x=5, y=0, z=0)), Atom(P(x=5, y=1, z=0)), Atom(P(x=5, y=2, z=0)), Atom(P(x=5, y=2, z=1)), Atom(P(x=6, y=0, z=0)), Atom(P(x=6, y=1, z=0)), Atom(P(x=6, y=2, z=0))]
================================================================================

Now we know that the lower bound on boxes is 31. But what does such a solution look like? Let us plot them using matplotlib, putting them all in a three-dimensional array

import matplotlib.pyplot as plt
import numpy as np

# Prepare some coordinates
x, y, z = np.indices((depth, depth, depth))

cubes = [((atom.object.x == x) & (atom.object.y == y) & (atom.object.z == z) )  for atom in atoms]

# Combine the objects into a single boolean array
voxelarray = cubes[0]
for cube in cubes[1:]:
  voxelarray |= cube

# Plot
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
ax.voxels(voxelarray, edgecolor='k')

ax.set(xticklabels=[], yticklabels=[], zticklabels=[])

plt.show()

And here is our trailer in glorius 3D:

Obviously the solution is not unique, as some stacks can be swapped and the whole configuration can be mirrored.

And just out of curiosity, if there were no gravity...

Solution with 26 boxes found
[Atom(P(x=0, y=0, z=0)), Atom(P(x=0, y=1, z=1)), Atom(P(x=0, y=2, z=1)), Atom(P(x=0, y=1, z=2)), Atom(P(x=0, y=2, z=2)), Atom(P(x=1, y=1, z=0)), Atom(P(x=1, y=2, z=0)), Atom(P(x=1, y=0, z=1)), Atom(P(x=1, y=0, z=2)), Atom(P(x=2, y=1, z=0)), Atom(P(x=2, y=2, z=0)), Atom(P(x=2, y=0, z=1)), Atom(P(x=2, y=0, z=2)), Atom(P(x=3, y=1, z=0)), Atom(P(x=3, y=2, z=0)), Atom(P(x=3, y=0, z=1)), Atom(P(x=3, y=0, z=2)), Atom(P(x=4, y=1, z=0)), Atom(P(x=4, y=2, z=0)), Atom(P(x=4, y=0, z=1)), Atom(P(x=5, y=1, z=0)), Atom(P(x=5, y=2, z=0)), Atom(P(x=5, y=0, z=1)), Atom(P(x=6, y=0, z=0)), Atom(P(x=6, y=1, z=0)), Atom(P(x=6, y=2, z=0))]
================================================================================
Solution with 25 boxes found
[Atom(P(x=0, y=2, z=0)), Atom(P(x=0, y=2, z=1)), Atom(P(x=0, y=0, z=2)), Atom(P(x=0, y=1, z=2)), Atom(P(x=1, y=2, z=0)), Atom(P(x=1, y=1, z=1)), Atom(P(x=1, y=0, z=2)), Atom(P(x=1, y=2, z=2)), Atom(P(x=2, y=2, z=0)), Atom(P(x=2, y=1, z=1)), Atom(P(x=2, y=0, z=2)), Atom(P(x=2, y=2, z=2)), Atom(P(x=3, y=2, z=0)), Atom(P(x=3, y=2, z=1)), Atom(P(x=3, y=0, z=2)), Atom(P(x=3, y=1, z=2)), Atom(P(x=4, y=1, z=0)), Atom(P(x=4, y=2, z=0)), Atom(P(x=4, y=0, z=1)), Atom(P(x=5, y=1, z=0)), Atom(P(x=5, y=2, z=0)), Atom(P(x=5, y=0, z=1)), Atom(P(x=6, y=0, z=0)), Atom(P(x=6, y=1, z=0)), Atom(P(x=6, y=2, z=0))]
================================================================================
Solution with 24 boxes found
[Atom(P(x=0, y=2, z=0)), Atom(P(x=0, y=2, z=1)), Atom(P(x=0, y=0, z=2)), Atom(P(x=0, y=1, z=2)), Atom(P(x=1, y=2, z=0)), Atom(P(x=1, y=1, z=1)), Atom(P(x=1, y=0, z=2)), Atom(P(x=1, y=2, z=2)), Atom(P(x=2, y=2, z=0)), Atom(P(x=2, y=1, z=1)), Atom(P(x=2, y=0, z=2)), Atom(P(x=2, y=2, z=2)), Atom(P(x=3, y=2, z=0)), Atom(P(x=3, y=1, z=1)), Atom(P(x=3, y=0, z=2)), Atom(P(x=4, y=1, z=0)), Atom(P(x=4, y=2, z=0)), Atom(P(x=4, y=0, z=1)), Atom(P(x=5, y=1, z=0)), Atom(P(x=5, y=2, z=0)), Atom(P(x=5, y=0, z=1)), Atom(P(x=6, y=0, z=0)), Atom(P(x=6, y=1, z=0)), Atom(P(x=6, y=2, z=0))]
================================================================================
Solution with 23 boxes found
[Atom(P(x=0, y=2, z=0)), Atom(P(x=0, y=2, z=1)), Atom(P(x=0, y=0, z=2)), Atom(P(x=0, y=1, z=2)), Atom(P(x=1, y=2, z=0)), Atom(P(x=1, y=1, z=1)), Atom(P(x=1, y=0, z=2)), Atom(P(x=1, y=2, z=2)), Atom(P(x=2, y=2, z=0)), Atom(P(x=2, y=1, z=1)), Atom(P(x=2, y=0, z=2)), Atom(P(x=3, y=2, z=0)), Atom(P(x=3, y=1, z=1)), Atom(P(x=3, y=0, z=2)), Atom(P(x=4, y=1, z=0)), Atom(P(x=4, y=2, z=0)), Atom(P(x=4, y=0, z=1)), Atom(P(x=5, y=1, z=0)), Atom(P(x=5, y=2, z=0)), Atom(P(x=5, y=0, z=1)), Atom(P(x=6, y=0, z=0)), Atom(P(x=6, y=1, z=0)), Atom(P(x=6, y=2, z=0))]
================================================================================
Solution with 22 boxes found
[Atom(P(x=0, y=2, z=0)), Atom(P(x=0, y=1, z=1)), Atom(P(x=0, y=0, z=2)), Atom(P(x=1, y=1, z=0)), Atom(P(x=1, y=0, z=1)), Atom(P(x=1, y=2, z=2)), Atom(P(x=2, y=0, z=0)), Atom(P(x=2, y=2, z=0)), Atom(P(x=2, y=0, z=1)), Atom(P(x=2, y=1, z=2)), Atom(P(x=3, y=0, z=0)), Atom(P(x=3, y=2, z=1)), Atom(P(x=3, y=1, z=2)), Atom(P(x=4, y=0, z=0)), Atom(P(x=4, y=2, z=0)), Atom(P(x=4, y=1, z=1)), Atom(P(x=5, y=1, z=0)), Atom(P(x=5, y=2, z=0)), Atom(P(x=5, y=0, z=1)), Atom(P(x=6, y=0, z=0)), Atom(P(x=6, y=1, z=0)), Atom(P(x=6, y=2, z=0))]
================================================================================
Solution with 21 boxes found
[Atom(P(x=0, y=1, z=0)), Atom(P(x=0, y=0, z=1)), Atom(P(x=0, y=2, z=2)), Atom(P(x=1, y=1, z=0)), Atom(P(x=1, y=2, z=1)), Atom(P(x=1, y=0, z=2)), Atom(P(x=2, y=0, z=0)), Atom(P(x=2, y=1, z=1)), Atom(P(x=2, y=2, z=2)), Atom(P(x=3, y=2, z=0)), Atom(P(x=3, y=0, z=1)), Atom(P(x=3, y=1, z=2)), Atom(P(x=4, y=0, z=0)), Atom(P(x=4, y=1, z=0)), Atom(P(x=4, y=2, z=1)), Atom(P(x=5, y=0, z=0)), Atom(P(x=5, y=2, z=0)), Atom(P(x=5, y=1, z=1)), Atom(P(x=6, y=0, z=0)), Atom(P(x=6, y=1, z=0)), Atom(P(x=6, y=2, z=0))]
================================================================================

Published: Fri 07 September 2018
Prolog Peano Puzzle

Finding Trade and Exchange Possibilities with Prolog

A friend recently had an interesting problem for me: There's a game for mobile phones where you have to exchange or transmute wares into each other, in order to create some strange supply chain.

Example of an earlier level.

A very complex later example.

The game is bleentoro, image copyright by them.

In this example (an earlier level) the goal is to create fish. We can transmute or exchange (I will only write about exchange from now on) an apple and a pair of cherries for a fish. Apple and cherry we have, or are created from thin air, in the context of the game. As you can see in the other example, this can get very complex.

The question is, given some input we have, in what order do we have to use the stations for exchange to create some output. How do we best model this problem?
Modeling the Problem

From a previous article we already know that we can encode numbers in Prolog terms using Peano numbers. Using these we can match terms that represent at least or exactly some number. Exactly 3 would be s(s(s(z))) while at least 3 would be s(s(s(X))) with the variable X representing the rest of the term for some $n$ so that $3 + x = n$.

Let's start with a small example: we have an aunt that gives us apples. We have bob who is willing to exchange an apple for two berries. And we have a cousin who will bake a cake if we give them two apples and five berries.

Let's formalize our rules in some way:
```
a: / -> 1A
b: A -> 2B
c: 2A, 5B -> 1C
```
We model this as some kind of inventory(A, B, C) with some slots for Peano encoded numbers. Meaning that inventory(s(z), z, s(C)) means an inventory of exactly one apple, zero berries, and at least one cake. Lucky us.

Let's define our possible exchanges as Prolog facts with exchange(Kind, Pre, Post) meaning we use the exchange method Kind and we get from inventory state Pre to inventory state Post. If we lose an item, we have it in Pre but don't in Post. If we gain one, it's the other way around. Of course, the rest of our inventory should stay the same.
```
exchange(a, inventory(A, B, C), inventory(s(A), B, C)).
exchange(b, inventory(s(A),B,C), inventory(A, s(s(B)), C)).
exchange(c, inventory(s(s(A)), s(s(s(s(s(B))))), C), inventory(A, B, s(C))).
```
As we can see, we can encode basically any numeric condition that uses only "just exactly" and "at least", as well as addition and subtraction and "setting" fields to specific values. If we have to give at least three apples, but all of them, to get a cake, the rule would be exchange(c, inventory(s(s(s(_))), B, C), inventory(z, B, s(C)))..
Programming the Solver

We want to create a predicate trades(Start, Goal, Path) where we give a Start inventory and a Goal inventory and want to get a Path of exchanges. The most obvious approach is the following:
```
trades(Goal, Goal, []).
trades(Start, Goal, [H|T]) :- exchange(H, Start, Temp), trades(Temp, Goal, T).
```
So we're done when our Goal is reached. And if not, we have to try some exchange. We try it and we quickly learn that...
```
?- trades(inventory(z, z, z), inventory(_, _, s(z)), Path).
ERROR: Out of local stack
```
... it doesn't work.

The problem is, that calling the exchange predicate in the recursive call always tries the first rule first and since that basically always works, we get into an infinite loop. Let's reorder the rules and put the generating rule first.
```
?- trades(inventory(z, z, z), inventory(_, _, s(z)), Path).
ERROR: Out of local stack
```
It still doesn't work. Let's look at the trace.
```
[trace]  ?- trades(inventory(z, z, z), inventory(_, _, s(z)), Path).
 Call: (8) trades(inventory(z, z, z), inventory(_700, _702, s(z)), _712) ? creep
 Call: (9) exchange(_960, inventory(z, z, z), _984) ? creep
 Exit: (9) exchange(a, inventory(z, z, z), inventory(s(z), z, z)) ? creep
 Call: (9) trades(inventory(s(z), z, z), inventory(_700, _702, s(z)), _962) ? creep
 Call: (10) exchange(_978, inventory(s(z), z, z), _1002) ? creep
 Exit: (10) exchange(b, inventory(s(z), z, z), inventory(z, s(s(z)), z)) ? creep
 Call: (10) trades(inventory(z, s(s(z)), z), inventory(_700, _702, s(z)), _980) ? creep
 Call: (11) exchange(_1000, inventory(z, s(s(z)), z), _1024) ? creep
 Exit: (11) exchange(a, inventory(z, s(s(z)), z), inventory(s(z), s(s(z)), z)) ? creep
 Call: (11) trades(inventory(s(z), s(s(z)), z), inventory(_700, _702, s(z)), _1002) ? creep
 Call: (12) exchange(_1018, inventory(s(z), s(s(z)), z), _1042) ? creep
 Exit: (12) exchange(b, inventory(s(z), s(s(z)), z), inventory(z, s(s(s(s(z)))), z)) ? creep
 Call: (12) trades(inventory(z, s(s(s(s(z)))), z), inventory(_700, _702, s(z)), _1020) ? creep
 Call: (13) exchange(_1040, inventory(z, s(s(s(s(z)))), z), _1064) ? creep
 Exit: (13) exchange(a, inventory(z, s(s(s(s(z)))), z), inventory(s(z), s(s(s(s(z)))), z)) ? creep
 Call: (13) trades(inventory(s(z), s(s(s(s(z)))), z), inventory(_700, _702, s(z)), _1042) ? creep
 Call: (14) exchange(_1058, inventory(s(z), s(s(s(s(z)))), z), _1082) ? creep
 Exit: (14) exchange(b, inventory(s(z), s(s(s(s(z)))), z), inventory(z, s(s(s(s(s(s(z)))))), z)) ? creep
 Call: (14) trades(inventory(z, s(s(s(s(s(s(z)))))), z), inventory(_700, _702, s(z)), _1060) ? creep
 Call: (15) exchange(_1080, inventory(z, s(s(s(s(s(s(z)))))), z), _1104) ? creep
 Exit: (15) exchange(a, inventory(z, s(s(s(s(s(s(z)))))), z), inventory(s(z), s(s(s(s(s(s(z)))))), z)) ? creep
 Call: (15) trades(inventory(s(z), s(s(s(s(s(s(z)))))), z), inventory(_700, _702, s(z)), _1082) ? creep
 Call: (16) exchange(_1098, inventory(s(z), s(s(s(s(s(s(z)))))), z), _1122) ? creep
```
We're still running into an infinite loop since we're just getting apples and trading them for berries.

There's a trick if our solution is a list. If the length(List, Length) predicate is called with both predicates free (length(L, _) for example), we get lists of increasing length as solutions. So if we can bind our solution to lists of increasing length, we won't get infinite loops of apple-berry-exchanges without ever trying short solutions.

So our query now looks like this:
```
?- length(Path,_), trades(inventory(z, z, z), inventory(_, _, s(z)), Path).
Path = [a, b, a, b, a, b, a, a, c] ;
Path = [a, b, a, b, a, a, b, a, c] ;
Path = [a, b, a, b, a, a, a, b, c] ;
Path = [a, b, a, a, b, b, a, a, c] ;
Path = [a, b, a, a, b, a, b, a, c] ;
Path = [a, b, a, a, b, a, a, b, c] ;
Path = [a, b, a, a, a, b, b, a, c] ;
Path = [a, b, a, a, a, b, a, b, c] ;
Path = [a, b, a, a, a, a, b, b, c] ...
```
As it turns out, there are a few ways to do it, but nothing shorter than 9 trades. We can even look, whether we have anything left over after we got our cake:
```
?- length(Path,_), trades(inventory(z, z, z), inventory(A, B, s(z)), Path).
Path = [a, b, a, b, a, b, a, a, c],
A = z,
B = s(z) .
```
Yeah. As it turns out, we have a berry to spare.
Conclusion

This is how you do it. And of course, you could simulate a shop with that, as one inventory slot could be the number of coins. If the inventory gets too large, you can split it into item types so you don't have to copy every unchanged item from the pre-condition to the post-condition. Then you can copy a term that might represent any number of items.

But beware: this algorithm (with Prolog's evaluation scheme) has not so good complexity. For a solution of length $n$ and $m$ rules it has a run time of $O(m^n)$. We create every possible solution list permutation but break early, if no rules are applicable.

Depending on our instantiation of Start and Goal, this algorithm doesn't terminate (if the system isn't terminating, the algorithm isn't either). Our system doesn't terminate as we can always apply a again and again. We're creating increasingly longer partial solutions which might turn out to be no solution at all. If there is no solution, you won't know. The algorithm just tries longer and longer lists trying to find a solution.

If you have any other conditions on your solution (some exchanges might be limited or necessary) you can put most of them after the trades call. Some conditions, like membership, can be put before.

Let's say we have a terminating system by taking away the a rule. If our Start is fully instantiated, then our algorithm terminates. Let's try another query and find out, with how many apples we need to start to get a cake:
```
?- trades(inventory(A, z, z), inventory(_, _, s(z)), Path).
ERROR: Out of local stack
```
Through the unification rules, it is possible that a rule creates terms on the left side of the rule. In this case, the algorithm tries to give us increasing amounts of berries to match that with the starting condition (only rule b is applied in reverse, always leaving us with one cake, which doesn't match the starting condition). There we need the length predicate.
```
?- length(Path,_), trades(inventory(A, z, z), inventory(_, _, s(z)), Path).
Path = [b, b, b, c],
A = s(s(s(s(s(_2334))))) .
```
We need at least five apples for a cake. But we knew that.

To recap:
- We saw again how useful it is, to encode the natural numbers in a term structure.
- A solver is easily written but the algorithm isn't fast.
- The rules can be applied in both directions because of Prolog's evaluation scheme. This leaves queries that would intuitively work divergent.
- By using the length predicate we can simulate breadth-first evaluation, while the default evaluation scheme is depth-first.
I think, there's a different approach using linear programming. I might try that in the future.
Published: Mon 27 August 2018
countdown haskell puzzle

Haskell Countdown Solver

I love watching 8 Out Of 10 Cats Does Countdown. It's a british panel show where a few comedians solve word puzzles and mathematics puzzles. It's based on the game show Countdown that's been running since 1982.

And I also like the host Jimmy Carr very much.

The Puzzles

The show consists of three different puzzles:

In the Letters round a team (or single contestant) picks 9 consonants and vowels from their respective piles. After the team chooses from one of both piles, the letter is revealed and put up on the board. They then choose again until 9 letters are on the board. After all letters are on the board, the teams have 30 seconds to form a word out of as many of the letters as possible. The team with the longest word scores that many points.

Letters Round

In the Conundrum a list of letters that already spell out a funny group of words is given by the host and the contestants try to find a 9 letter anagram. The first team to buzzer in the correct answer gains 10 points. This is the final round of the game.

Conundrum

In the Numbers round a team (or single contestant) picks 6 numbers from two piles. A random number is picked and the contestants have to create a mathematical formula out of the picked numbers that equal the target number. The specific rules follow.

Numbers Round
Numbers Round

According to Wikipedia there are 6 columns of 4 numbers (4 rows). One column contains the numbers 25, 50, 75, and 100 and the other 5 columns contain the numbers 1 to 10 twice each. That's not exactly the same as on the show, as there appear to be 7 columns of 4 numbers but not every stack has all 4 number cards.

Contestants choose a total of 6 numbers from the "small" and "large" columns. The usual pick would be 2 to 3 "large" numbers and the rest small ones.

Then a random number between 100 und 999 is generated. The contestants try to create a mathematical formula that has the random number as the result. In the formulare only the following operations are allowed:
- Addition
- Multiplication
- Subtraction (that does not result in a negative number)
- Division (that results in an integer)
The contestants do not have to use all the numbers to get the target number.

Finished Numbers Round

Since I am really not good at this, I wrote a solver for the numbers game in Haskell.
Solver

How do we go about this problem? First of all, we have to decide on a data structure. Let's start with a traditional term structure where a Term is a binary operation and two terms, or an integer value. An Operation is the function mapping two integers to a third, and a string to pretty-print the operation. The mapping is not total and might fail.
```
data Term = Op Operation Term Term | Val Integer
data Operation = Operation {fun :: (Integer -> Integer -> Maybe Integer), name :: String}

instance Show Term where
  show (Val x) = show x
  show (Op op l r) = "(" ++ (show l) ++ (name op) ++ show r ++ ")"
```
Now we have to define our operations and what they should do. Since we're using Maybe, and it's a monad, we can use monadic notation where it suits us:
```
operations = [  Operation (\x y -> return (x + y)) "+"
             ,  Operation (\x y -> return (x * y)) "*"
             ,  Operation (\x y -> guard (x > y) >> return (x - y)) "-"
             ,  Operation (\x y -> guard (x `mod` y == 0) >> return (x `div` y)) "/"
             ]
eval :: Term -> Maybe Integer
eval (Val v) = Just v
eval (Op o l r) = do l' <- eval l
                     r' <- eval r
                     (fun o) l' r'
```
As we can see, the preconditions to the subtraction and division are encoded using the guard function. If the function returns true, the whole operation returns Nothing, otherwise it returns Just the value.

Since Haskell uses lazy evaluation, we can define a list of all possible terms given a list of integers and filter for the ones that evaluate to the correct number.

As we create a Term we have to be careful to not use any of the numbers again. So we will split our numbers in two sets. One that can be used in the left Term of an Operation and the rest can be used in the right Term. We will create all terms using our definition. A Term is either the Value of one of the numbers for that subterm, or an operation with some numbers reserved for the left term and some numbers reserved for the right term.
```
terms :: [Integer] -> [Term]
terms nums =  do  n <- nums -- choose a number
                  [Val n] -- one solution
              ++
              do  op <- operations -- choose an operation
                  (l, r) <- subset_split nums -- choose a split
                  guard $ l /= []
                  guard $ r /= []
                  ls <- (terms l) -- choose a term for the left side
                  guard (eval ls /= Nothing)
                  rs <- (terms r) -- choose a term for the right side
                  guard (eval rs /= Nothing)
                  [Op op ls rs] -- one solution
```
For performance reasons we already evaluate the generated subterms in order to check whether the operations are valid. There are a lot of terms that can not be evaluated and we want to throw them out as early as possible (or not even generate them).

We haven't defined subset_split yet. The idea is, to just split the list of numbers in two groups in every possible way:
```
subsets :: [Integer] -> [[Integer]]
subsets []  = [[]]
subsets (x:xs) = subsets xs ++ map (x:) (subsets xs)

subset_split nums = do l <- subsets nums
                       return (l, nums \\ l)
```
Now all that remains is to generate all terms given a set of numbers, evaluate them, filter out the ones that don't have the correct value, and return the first one:
```
solve :: [Integer] -> Integer -> Maybe Term
solve nums target = let possible_solutions = map (\x -> (x, eval x)) (terms nums)
                        solutions = filter (\x -> (snd x) == (Just target)) possible_solutions
                    in listToMaybe $ map fst solutions
```
And in our main we only read the list of numbers and the target as arguments and print the solution:
```
main = do args <- getArgs
          let nums = (read (args !! 0))::[Integer]
              target = (read (args !! 1))::Integer
          print $ solve nums target
```
And there we go:
```
$ ./countdown "[100, 75, 2, 10, 3, 8]" 746
Just (8+(3+((75*(100-2))/10)))
```
Our solution is $8+3+{{75*(100-2)} \over 10}$ while the solution from the show was $10*75 - {8 \over 2}$.

As you can see, we do not always get the easiest solution but we usually have a solution quite fast. This one took about a tenth of a second. The 30 seconds given to the contestant are usually quite enough. To make sure we had the "nicest" solution, we probably would need to generate all terms and sort them by size. The smallest one should be the nicest one.
Published: Mon 21 March 2016
python backtracking puzzle

Solving Tatham's Puzzles - Signpost (backtracking)

I've started writing solvers for games of Simon Tatham's Portable Puzzle Collection. Those Problems can usually be reduced to some NP-complete problem. This is quite a fun exercise.

Unsolved Signpost Puzzle

We will be continuing with the puzzle "Signpost", which is a variant of Hamiltonian paths in a directed graph where, for some fields or nodes, the position within the path is given.

The goal of the puzzle is to connect the Node "1" to some node, which is then node "2" in the direction the arrow points to. All numbers need to be connected to their successor until all nodes have a single unique number from "1" to (in this example) "25".

We will be solving this problem in two different ways, both with the help of the programming language Python. But first of all we'll need to convert a puzzle into a machine readable format. Luckily, Signpost has an easy and accessible format. The problem instance from the image is the instance 5x5:1cceefcfggeeccghcac3e12hch10ah25a. It starts off with the size of the problem and then a list of letters that show the directions of the arrow where a is north or 12 o'clock and all other letters go clockwise with h being north-west. A letter can be prefixed by a number which is the hint for that node. That node's number is fixed. The cells or nodes are defined line by line from the top left to the bottom right.

Problems don't need to start with 1 or end with the last number. The start and the end of the path can be anywhere but puzzles generated with the ends at opposite corner look nicer.

First we define what our data looks like. We implement a simple form of directed graph that we define as a set of nodes and their successors. From that we derive following definition:
```
class Node(object):
    hint = None
    following = frozenset()

    def __init__(self, direction):
        self.direction = direction
```
We don't want to use any of that code in a library or something something so we don't do anything fancy here and we use the simplest class definition that comes to mind.

^(?P<size_x>\d+)x(?P<size_y>\d+)$

(?P<defn>\d*[a-h])

Now we go on to parsing the problem. The parsing function should get the problem as a string and return a 3-tuple containing width and height of the problem as well as a mapping of positions within the grid to ``Node`` instances.

Onto actually parsing. We can split at the colon and then regular expressions to the rescue. The size definition has the form ^(?P<size_x>\d+)x(?P<size_y>\d+)$ while a single node definition has the form (?P<defn>\d*[a-h]). As we can see, the digits are optional but we need to greedily accept them. All non-overlapping occurences in order will be our cell definitions.

In code we, again, do nothing too too fancy:
```
def parse(puzzle_str):
    directions = {
        'a': (0,-1), 'b': (1,-1), 'c': (1,0), 'd': (1,1),
        'e': (0,1),  'f': (-1,1), 'g': (-1,0),'h': (-1,-1)
    }

    size, _, definition = puzzle_str.partition(':')
    r_size = re.compile(r"^(?P<size_x>\d+)x(?P<size_y>\d+)$")
    r_defn = re.compile(r"(?P<defn>\d*[a-h])")
    size_t = tuple(map(int, r_size.match(size).groups()))
    w, h = size_t
    nodes = {}
    for n, m in enumerate(r_defn.finditer(definition)):
        pos = (n % w, n // w)
        nodes[pos] = Node(directions[m.group(0)[-1:]])
        hint = m.group(0)[:-1]
        if hint:
            nodes[pos].hint = int(hint)
    return w, h, nodes
```
And we also need to fill in the successor Nodes. This is no problem at all.
```
w, h, nodes = parse(sys.argv[1])
from itertools import product
for x, y in product(range(w), range(h)):
    this_node = nodes[(x,y)]
    dx, dy = this_node.direction
    x, y = x + dx, y + dy
    while x >= 0 and x < w and y >= 0 and y < h:
        this_node.following = this_node.following | frozenset([nodes[x,y]])
        x, y = x + dx, y + dy
```
And now on to actually solving this thing. The idea is, that we slowly build a list of Nodes that is our solution. The first element is the Node with the hint "1" and we iterate through all nodes that follow that Node and do not have a hint that is other than "2", and so on. And of course, the Node that we add to the solution is not allowed to have already been part of the solution and if there is a node with that hint, it should be this exact one. If we find such a conflict in our solution, we reject that solution and return to an earlier partial solution.

Instead of using recursion, we build a queue of partial solutions until we get a complete solution. This is functionally identical to restarting a solver function with partial solutions of increasing size. In that case, the call stack manages the backtracking. But we'll do it with a queue this time. We build the queue so that every element in the list is either
1. a partial solution (beginning from 1) without conflicts
2. a partial solution (beginning from 1) where the last element is in some conflict
3. a complete solution
For case 1, we remove the partial solution that has the length n and add all partial solutions of length n+1 by iterating through all successor Nodes that have not yet been part of the partial solution. For case 2, we reject that partial solution. For case 3, we immediately return that solution and are finished.
```
def solve_dhc(nodes):
    # get the Node that has hint 1
    first = next(iter(n for n in nodes.values() if n.hint == 1))
    # get the set of all hints, so it is easier to check for conflicts
    hints = frozenset(n.hint for n in nodes.values() if n.hint is not None)
    from collections import deque
    q = deque([[first]]) # initializing queue
    while q:
        curr = q.pop()
        idx = len(curr)
        if (idx in hints or curr[-1].hint is not None) and curr[-1].hint != idx:
            continue # case 2
        if idx == len(nodes):
            return curr # case 3
        for _next in curr[-1].following: # case 1
            if _next not in curr:
                q.append(curr + [_next])
```
This algorithm terminates because we always remove one element from the queue of some length n and possibly add a large but finite amount of elements of length n+1 and the algorithm terminates when we have an element in the queue of some possibly large but finite length.

Sidenote: if you use popleft instead of pop you get a breadth-first search instead the depth-first search, that is implemented here. Add to the right and pop from the right for depth-first and add to the right and pop from the left for breadth-first. Since every proper solution has the exact same length/depth (that we know and never go further anyways), searching depth-first is strictly better than breadth-first.

And given a solution, we still have to print it. That's the easy part.
```
dhc = solve_dhc(nodes)

fill = len(str(w*h))
for l in range(h):
    for c in range(w):
        node = nodes[(c,l)]
        print(str(dhc.index(node) + 1).rjust(fill), end=" ")
    print("")
```
Solved Puzzle

And there's the solution for our puzzle.
```
$ python3 signpost_dhc.py 5x5:1cceefcfggeeccghcac3e12hch10ah25a
 1 20  9  2 21
23 14 13 22 24
15  5  7  6  8
18 19 11  3 12
16 17 10  4 25
```