A case study in concise code - Algorithm X

Suppose we have a set U = {1, 2, 3, 4, 5, 6, 7} and a set S of certain subsets of U:

S =

{

    {1, 4, 7},

    {1, 4},

    {4, 5, 7},

    {3, 5, 6},

    {2, 3, 6, 7},

    {2, 7}

}

Each element of S is some subset of U.

Suppose we start taking some arbitrary subsets of S, like this:

{

    {1, 4, 7},

    {1, 4},

}

{

    {1, 4},

    {3, 5, 6},

    {2, 3, 6, 7},

}

...

How can we find a subset of S we will call S* such that every element of U appears in exactly one element of S*.

This is an example of the exact cover problem.

In our case the solution is this:

S* =

{

    {1, 4},

    {3, 5, 6},

    {2, 7}

}

Every element of U = {1, 2, 3, 4, 5, 6, 7} appears in exactly one element of S*.

An example of a non-solution is:

{

    {2, 3, 6, 7},

}

Not all elements of U appear in an element of the set.

Another example of a non-solution is:

{

    {1, 4, 7},

    {3, 5, 6},

    {2, 3, 6, 7},

    {2, 7}

}

Some elements of U appear in more than one element of the set.

---

Simple Algorithm X

Algorithm X is an algorithm developed by Donald Knuth that solves the exact cover problem. Let's implement Algorithm X in Haskell.

Now, don't panic. You don't necessarily need to know Haskell, as I will describe everything that's happening.

We can later use this to solve Sudoku puzzles (in another blog post).

Here come some imports.

module AlgorithmX where

import Data.IntSet hiding (map, filter, foldl)
import Data.IntMap (IntMap)
import qualified Data.IntMap as IntMap
import Data.List

import Text.Tabular hiding (row)
import Text.Tabular.AsciiArt

Don't panic.

Let's get some sets going.

set0 :: IntSet
set0 = fromList [1, 4, 7]

set0 is an IntSet (a set of Ints). Internally, it is represented as an immutable tree.

We can interact with it in the ghci interpreter:

ghci> set0
fromList [1,4,7]

ghci> size set0
3

ghci> :t member
member :: Key -> IntSet -> Bool

The function member takes a Key and an IntSet and returns a Bool indicating whether the keys is a member of the set.

ghci> :i Key
type Key :: *
type Key = Int
        -- Defined in `Data.IntSet.Internal'

Key is a type synonym for Int.

Checking element membership:

ghci> set0
fromList [1,4,7]

ghci> member 1 set0
True

ghci> member 0 set0
False

Let's get some more sets going:

sets1 :: [IntSet]
sets1 = map fromList [ [1, 4, 7]
                     , [1, 4]
                     , [4, 5, 7]
                     , [3, 5, 6]
                     , [2, 3, 6, 7]
                     , [2 ,7]
                     ]

sets1 is a list of IntSets.

sets1 is defined as a map that applies the function fromList to every element of the given list. (I know, Haskellers like to format lists like that... don't worry about it.)

I've prepared some printing functions earlier so we can print our data structures.

• printList
• printScan
• instance Show SparseMatrix

(Definitions in Appendix.)

Here's sets1:

ghci> printList "Set" sets1
+---++--------------------+
|   ||                Set |
+===++====================+
| 0 ||   fromList [1,4,7] |
| 1 ||     fromList [1,4] |
| 2 ||   fromList [4,5,7] |
| 3 ||   fromList [3,5,6] |
| 4 || fromList [2,3,6,7] |
| 5 ||     fromList [2,7] |
+---++--------------------+

In Algorithm X, we think of our sets as a matrix, like this:

	1	2	3	4	5	6	7
0	1	0	0	1	0	0	1
1	1	0	0	1	0	0	0
2	0	0	0	1	1	0	1
3	0	0	1	0	1	1	0
4	0	1	1	0	0	1	1
5	0	1	0	0	0	0	1

The column labels are elements of U.
Each row represents an elements of S, having the value 1 denoting it contains that element of U or 0 denoting it does not contain it. So row 0 is the set [1, 4, 7]: you can see 1s in columns 1, 4 and 7 and 0s in the other columns.

In order to achieve exact cover we need to pick out the collection of rows such that exactly one 1 appears in every column of the matrix. We can use the matrix to visualise which sets overlap by looking at the columns. Additionally, if we search for a solution by crossing out rows we don't want to use, we can identify a dead end by finding a column with only 0s.

Now, this matrix will be sparse. Let's make a sparse matrix. Our sparse matrix will have rows but we can keep the IntSet representation of each row.

type Row = IntSet

Row is type synonym for IntSet.

ghci> set0 :: IntSet
fromList [1,4,7]

ghci> set0 :: Row   
fromList [1,4,7]

Our interpreter accepts that set0 is an IntSet and also a Row.

We would like to have a collection of these rows.

type Rows = IntMap Row

Rows is an associative mapping from an Int to a Row. IntMap is internally represented as an immutable tree.

In order to make Rows, we need to give our Rows some indices.

ghci> printList "(Key, Row)" $ zip [0..] sets1
+---++------------------------+
|   ||             (Key, Row) |
+===++========================+
| 0 ||   (0,fromList [1,4,7]) |
| 1 ||     (1,fromList [1,4]) |
| 2 ||   (2,fromList [4,5,7]) |
| 3 ||   (3,fromList [3,5,6]) |
| 4 || (4,fromList [2,3,6,7]) |
| 5 ||     (5,fromList [2,7]) |
+---++------------------------+

Now we can put the rows into an IntMap of Rows.

rows1 = IntMap.fromList $ zip [0..] sets1

ghci> printList "(Key, Row)" $ IntMap.toList rows1                  
+---++------------------------+
|   ||             (Key, Row) |
+===++========================+
| 0 ||   (0,fromList [1,4,7]) |
| 1 ||     (1,fromList [1,4]) |
| 2 ||   (2,fromList [4,5,7]) |
| 3 ||   (3,fromList [3,5,6]) |
| 4 || (4,fromList [2,3,6,7]) |
| 5 ||     (5,fromList [2,7]) |
+---++------------------------+

ghci> rows1 ! 0
fromList [1,4,7]

ghci> rows1 ! 3
fromList [3,5,6]

We can query rows by their key (in O(log n)).

We can remove a row:

ghci> printList "(Key, Row)" $ IntMap.toList $ IntMap.delete 3 rows1
+---++------------------------+
|   ||             (Key, Row) |
+===++========================+
| 0 ||   (0,fromList [1,4,7]) |
| 1 ||     (1,fromList [1,4]) |
| 2 ||   (2,fromList [4,5,7]) |
| 3 || (4,fromList [2,3,6,7]) |
| 4 ||     (5,fromList [2,7]) |
+---++------------------------+

Without affecting the original rows1:

ghci> printList "(Key, Row)" $ IntMap.toList rows1
+---++------------------------+
|   ||             (Key, Row) |
+===++========================+
| 0 ||   (0,fromList [1,4,7]) |
| 1 ||     (1,fromList [1,4]) |
| 2 ||   (2,fromList [4,5,7]) |
| 3 ||   (3,fromList [3,5,6]) |
| 4 || (4,fromList [2,3,6,7]) |
| 5 ||     (5,fromList [2,7]) |
+---++------------------------+

The Rows with the deleted row are a new data structure. Both the old and the new are immutable and the new re-uses much of the old structure under the hood.

This will be handy later when we will want to delete rows and columns from our sparse matrix and then backtrack.

We want to use something like Rows to represent our sparse matrix, but at the moment when we query for the existence of elements within our rows, we cannot tell whether the queried element is within the matrix or not.

ghci> rows1 ! 0
fromList [1,4,7]

ghci> member 1 $ rows1 ! 0
True

1 is a member of row 0.

ghci> member 2 $ rows1 ! 0
False

2 is a not a member of row 0.

ghci> member 10 $ rows1 ! 0
False

10 is not a member of row 0. But is it within the bounds of the matrix?

We need another piece of information to keep track of what is in our sparse matrix: which columns are in the matrix?

type ActiveCols = IntSet

ActiveCols is a type synonym for IntSet.

In our case, we would want to use the unions of our rows1 to form the active columns.

ghci> printList "(Key, Row)" $ IntMap.toList rows1
+---++------------------------+
|   ||             (Key, Row) |
+===++========================+
| 0 ||   (0,fromList [1,4,7]) |
| 1 ||     (1,fromList [1,4]) |
| 2 ||   (2,fromList [4,5,7]) |
| 3 ||   (3,fromList [3,5,6]) |
| 4 || (4,fromList [2,3,6,7]) |
| 5 ||     (5,fromList [2,7]) |
+---++------------------------+

ghci> unions rows1
fromList [1,2,3,4,5,6,7]

data SparseMatrix = SparseMatrix Rows ActiveCols

SparseMatrix is a new data type containing Rows and ActiveCols.

m1 = SparseMatrix rows1 (unions rows1)

ghci> print m1

+---++---+---+---+---+---+---+---+
|   || 1 | 2 | 3 | 4 | 5 | 6 | 7 |
+===++===+===+===+===+===+===+===+
| 0 || 1 | 0 | 0 | 1 | 0 | 0 | 1 |
| 1 || 1 | 0 | 0 | 1 | 0 | 0 | 0 |
| 2 || 0 | 0 | 0 | 1 | 1 | 0 | 1 |
| 3 || 0 | 0 | 1 | 0 | 1 | 1 | 0 |
| 4 || 0 | 1 | 1 | 0 | 0 | 1 | 1 |
| 5 || 0 | 1 | 0 | 0 | 0 | 0 | 1 |
+---++---+---+---+---+---+---+---+

The matrix prints out like this, but internally it is just Rows (i.e. IntMap IntSet) and ActiveCols (i.e. IntSet).

In order to search for an exact cover solution, we should pick a column and try to fill it. For the moment, column 1 seems as good as any. We have two choices of how to satisfy column 1: with row 0 or row 1. Row 0 seems as good as any.

Now we have to remove from the matrix all rows that conflict with row 0. Row 0 contains elements 1, 4 and 7. Any other rows that use these elements must be removed, otherwise we will have some elements of U in multiple elements of S*. In other words, we must only keep rows disjoint with row 0.

ghci> :set +m

(Set multiline mode.)

ghci> (SparseMatrix (IntMap.filter ((rows1 ! 0) `disjoint`) rows1)
ghci|               (unions rows1))

+---++---+---+---+---+---+---+---+
|   || 1 | 2 | 3 | 4 | 5 | 6 | 7 |
+===++===+===+===+===+===+===+===+
| 3 || 0 | 0 | 1 | 0 | 1 | 1 | 0 |
+---++---+---+---+---+---+---+---+

Only row 3 is disjoint with row 0.

Another thing we need to do is remove columns from the matrix that we satisfied by picking row 0. We remove these from the ActiveCols of the matrix.

ghci> unions rows1
fromList [1,2,3,4,5,6,7]

ghci> rows1 ! 0
fromList [1,4,7]

ghci> unions rows1 `difference` (rows1 ! 0)
fromList [2,3,5,6]

ghci> (SparseMatrix (IntMap.filter ((rows1 ! 0) `disjoint`) rows1)
ghci|               (unions rows1 `difference` (rows1 ! 0)))

+---++---+---+---+---+
|   || 2 | 3 | 5 | 6 |
+===++===+===+===+===+
| 3 || 0 | 1 | 1 | 1 |
+---++---+---+---+---+

Unfortunately, we can see that column 2 cannot be satisfied by this remaining row. So row 0 cannot be part of the solution. We should go back and pick another row.

Let's go back and try row 1.

ghci> print m1

+---++---+---+---+---+---+---+---+
|   || 1 | 2 | 3 | 4 | 5 | 6 | 7 |
+===++===+===+===+===+===+===+===+
| 0 || 1 | 0 | 0 | 1 | 0 | 0 | 1 |
| 1 || 1 | 0 | 0 | 1 | 0 | 0 | 0 |
| 2 || 0 | 0 | 0 | 1 | 1 | 0 | 1 |
| 3 || 0 | 0 | 1 | 0 | 1 | 1 | 0 |
| 4 || 0 | 1 | 1 | 0 | 0 | 1 | 1 |
| 5 || 0 | 1 | 0 | 0 | 0 | 0 | 1 |
+---++---+---+---+---+---+---+---+

ghci> (SparseMatrix (IntMap.filter ((rows1 ! 1) `disjoint`) rows1)
ghci|               (unions rows1 `difference` (rows1 ! 1)))

+---++---+---+---+---+---+
|   || 2 | 3 | 5 | 6 | 7 |
+===++===+===+===+===+===+
| 3 || 0 | 1 | 1 | 1 | 0 |
| 4 || 1 | 1 | 0 | 1 | 1 |
| 5 || 1 | 0 | 0 | 0 | 1 |
+---++---+---+---+---+---+

This matrix represents the remaining options after having picked row 1. (Again, other rows conflicting with the choice and the filled columns have been removed.)

No columns have all 0s, so we could keep searching within this matrix in hopes of finding the complete solution.

After playing around like this, we are ready to look at Algorithm X and understand what it is trying to do.

From Wikipedia:

• If the matrix A has no columns, the current partial solution is a valid solution; terminate successfully.
  • Otherwise choose a column c (deterministically).
  • Choose a row r such that Ar, c = 1 (nondeterministically).
  • Include row r in the partial solution.
  • For each column j such that Ar, j = 1,
    • for each row i such that Ai, j = 1,
      • delete row i from matrix A.
    • delete column j from matrix A.
  • Repeat this algorithm recursively on the reduced matrix A.

"Deterministically" here means we have a procedure of picking some column right away, and "nondeteministically" means that the choice involves some kind of space search.

We have, in fact, just performed one iteration of Algorithm X manually. We chose a column, chose a row and reduced the matrix.

Instead of performing these steps through iteration

• ...
  • For each column j such that Ar, j = 1,
    • for each row i such that Ai, j = 1,
      • delete row i from matrix A.

we used equivalent functionality to filter the rows by the criterion of being disjoint with the selected row (the disjoint operation runs in O(n+m)), and instead of

• ...
  • For each column j such that Ar, j = 1,
    • ...
    • delete column j from matrix A.

we used the difference between the active columns and the selected row (also O(n+m)). (Keep in mind our operations do not mutate the original data structure and allow easy backtracking.)

It's time to think about writing a recursive function.

-- NOT final version

type IsCompleteSolution = Bool

scanAlgoXSimple' :: SparseMatrix
                 -> [Key]
                 -> [([Key], SparseMatrix, IsCompleteSolution)]
...

scanAlgoXSimple' is a function taking a matrix and a partial solution (list of Keys) and returning a list of intermediate results: a list of state tuples of (a solution, a matrix, a boolean indicating if the solution is complete or partial).

-- NOT final version

type IsCompleteSolution = Bool

scanAlgoXSimple' :: SparseMatrix
                 -> [Key]
                 -> [([Key], SparseMatrix, IsCompleteSolution)]
scanAlgoXSimple' m@(SparseMatrix rows activeCols) solution
    | size activeCols == 0 = [(solution, m, True)]
    | otherwise            = (solution, m, False) : []

• If there are no more active columns, scanAlgoXSimple' returns the solution, the matrix, and True indicating that this a complete solution.
• Otherwise, for the moment we will just return the current partial solution, the current matrix and False indicating that this is a partial solution. And then we don't do anything further.

ghci> printScan $ scanAlgoXSimple' m1 []
([],
+---++---+---+---+---+---+---+---+
|   || 1 | 2 | 3 | 4 | 5 | 6 | 7 |
+===++===+===+===+===+===+===+===+
| 0 || 1 | 0 | 0 | 1 | 0 | 0 | 1 |
| 1 || 1 | 0 | 0 | 1 | 0 | 0 | 0 |
| 2 || 0 | 0 | 0 | 1 | 1 | 0 | 1 |
| 3 || 0 | 0 | 1 | 0 | 1 | 1 | 0 |
| 4 || 0 | 1 | 1 | 0 | 0 | 1 | 1 |
| 5 || 0 | 1 | 0 | 0 | 0 | 0 | 1 |
+---++---+---+---+---+---+---+---+
,False)

The solution is empty, the matrix is full, and the solution is partial. And we don't go further.

-- NOT final version

type IsCompleteSolution = Bool

scanAlgoXSimple' :: SparseMatrix
                 -> [Key]
                 -> [([Key], SparseMatrix, IsCompleteSolution)]
scanAlgoXSimple' m@(SparseMatrix rows activeCols) solution
    | size activeCols == 0 = [(solution, m, True)]
    | otherwise =
        let (r, row)  = IntMap.findMin rows
            m'        = SparseMatrix (IntMap.filter (row `disjoint`) rows)
                                     (activeCols `difference` row)
        in  (solution, m, False) : scanAlgoXSimple' m' (r:solution)

Here we try taking the first available row. Then we reduce the matrix as we did before: the rows are filtered for being disjoint with the selected row and the row's 1s are removed from the active columns with difference - giving us the reduced matrix m'. We return the current state tuple and recursively return the rest of the states by calling scanAlgoXSimple' m' (r:solution). We pass the reduced matrix to scanAlgoXSimple', as well as prepending the selected row key r to the current partial solution. (Our solution list will grow backwards.)

ghci> printScan $ scanAlgoXSimple' m1 []
([],
+---++---+---+---+---+---+---+---+
|   || 1 | 2 | 3 | 4 | 5 | 6 | 7 |
+===++===+===+===+===+===+===+===+
| 0 || 1 | 0 | 0 | 1 | 0 | 0 | 1 |
| 1 || 1 | 0 | 0 | 1 | 0 | 0 | 0 |
| 2 || 0 | 0 | 0 | 1 | 1 | 0 | 1 |
| 3 || 0 | 0 | 1 | 0 | 1 | 1 | 0 |
| 4 || 0 | 1 | 1 | 0 | 0 | 1 | 1 |
| 5 || 0 | 1 | 0 | 0 | 0 | 0 | 1 |
+---++---+---+---+---+---+---+---+
,False)

([0],
+---++---+---+---+---+
|   || 2 | 3 | 5 | 6 |
+===++===+===+===+===+
| 3 || 0 | 1 | 1 | 1 |
+---++---+---+---+---+
,False)

([3,0],
+--++---+
|  || 2 |
+==++===+
+--++---+
,False)*** Exception: findMin: empty map has no minimal element

We picked row 0. Next we only had one choice, to pick row 3. Next we were left with no rows but column 2 had not been satisfied. We have not handled the case where are are still active columns but no rows, so we got an error.

-- Final version

type IsCompleteSolution = Bool

scanAlgoXSimple' :: SparseMatrix
                 -> [Key]
                 -> [([Key], SparseMatrix, IsCompleteSolution)]
scanAlgoXSimple' m@(SparseMatrix rows activeCols) solution
    | size activeCols == 0 = [(solution, m, True)]
    | otherwise =
        (solution, m, False)
            : [ s | (r, row) <- IntMap.toList rows,
                    let m' = SparseMatrix (IntMap.filter (row `disjoint`) rows)
                                          (activeCols `difference` row),
                    s <- scanAlgoXSimple' m' (r:solution) ]

Instead of picking a row out of the matrix as before, we now iterate over all rows ((r, row) <- IntMap.toList rows). For each row we reduce the matrix and recursively return all state results for that selection (s <- scanAlgoXSimple' m' (r:solution)). When we go on to the next row, the matrix will be reduced according to that row selection and the results of s <- scanAlgoXSimple' m' (r:solution) with the new row will keep being put into the same list (it's all happening in one list comprehension), and so on.

ghci> printScan $ scanAlgoXSimple' m1 []
([],
+---++---+---+---+---+---+---+---+
|   || 1 | 2 | 3 | 4 | 5 | 6 | 7 |
+===++===+===+===+===+===+===+===+
| 0 || 1 | 0 | 0 | 1 | 0 | 0 | 1 |
| 1 || 1 | 0 | 0 | 1 | 0 | 0 | 0 |
| 2 || 0 | 0 | 0 | 1 | 1 | 0 | 1 |
| 3 || 0 | 0 | 1 | 0 | 1 | 1 | 0 |
| 4 || 0 | 1 | 1 | 0 | 0 | 1 | 1 |
| 5 || 0 | 1 | 0 | 0 | 0 | 0 | 1 |
+---++---+---+---+---+---+---+---+
,False)

([0],
+---++---+---+---+---+
|   || 2 | 3 | 5 | 6 |
+===++===+===+===+===+
| 3 || 0 | 1 | 1 | 1 |
+---++---+---+---+---+
,False)

([3,0],
+--++---+
|  || 2 |
+==++===+
+--++---+
,False)

([1],
+---++---+---+---+---+---+
|   || 2 | 3 | 5 | 6 | 7 |
+===++===+===+===+===+===+
| 3 || 0 | 1 | 1 | 1 | 0 |
| 4 || 1 | 1 | 0 | 1 | 1 |
| 5 || 1 | 0 | 0 | 0 | 1 |
+---++---+---+---+---+---+
,False)

([3,1],
+---++---+---+
|   || 2 | 7 |
+===++===+===+
| 5 || 1 | 1 |
+---++---+---+
,False)

([5,3,1],
+--++--+
|  ||  |
+==++==+
+--++--+
,True)

([4,1],
+--++---+
|  || 5 |
+==++===+
+--++---+
,False)

([5,1],
+---++---+---+---+
|   || 3 | 5 | 6 |
+===++===+===+===+
| 3 || 1 | 1 | 1 |
+---++---+---+---+
,False)

([3,5,1],
+--++--+
|  ||  |
+==++==+
+--++--+
,True)

([2],
+--++---+---+---+---+
|  || 1 | 2 | 3 | 6 |
+==++===+===+===+===+
+--++---+---+---+---+
,False)

([3],
+---++---+---+---+---+
|   || 1 | 2 | 4 | 7 |
+===++===+===+===+===+
| 0 || 1 | 0 | 1 | 1 |
| 1 || 1 | 0 | 1 | 0 |
| 5 || 0 | 1 | 0 | 1 |
+---++---+---+---+---+
,False)

([0,3],
+--++---+
|  || 2 |
+==++===+
+--++---+
,False)

([1,3],
+---++---+---+
|   || 2 | 7 |
+===++===+===+
| 5 || 1 | 1 |
+---++---+---+
,False)

([5,1,3],
+--++--+
|  ||  |
+==++==+
+--++--+
,True)

([5,3],
+---++---+---+
|   || 1 | 4 |
+===++===+===+
| 1 || 1 | 1 |
+---++---+---+
,False)

([1,5,3],
+--++--+
|  ||  |
+==++==+
+--++--+
,True)

([4],
+---++---+---+---+
|   || 1 | 4 | 5 |
+===++===+===+===+
| 1 || 1 | 1 | 0 |
+---++---+---+---+
,False)

([1,4],
+--++---+
|  || 5 |
+==++===+
+--++---+
,False)

([5],
+---++---+---+---+---+---+
|   || 1 | 3 | 4 | 5 | 6 |
+===++===+===+===+===+===+
| 1 || 1 | 0 | 1 | 0 | 0 |
| 3 || 0 | 1 | 0 | 1 | 1 |
+---++---+---+---+---+---+
,False)

([1,5],
+---++---+---+---+
|   || 3 | 5 | 6 |
+===++===+===+===+
| 3 || 1 | 1 | 1 |
+---++---+---+---+
,False)

([3,1,5],
+--++--+
|  ||  |
+==++==+
+--++--+
,True)

([3,5],
+---++---+---+
|   || 1 | 4 |
+===++===+===+
| 1 || 1 | 1 |
+---++---+---+
,False)

([1,3,5],
+--++--+
|  ||  |
+==++==+
+--++--+
,True)

We have found all of the solutions. They are all permutations of the same solution {1, 3, 5}.

The output is a little too long, so let's write a function to stop after the first solution is found.

upToComplete :: [([Key], SparseMatrix, IsCompleteSolution)]
             -> [([Key], SparseMatrix, IsCompleteSolution)]
upToComplete states =
    case break (\(_, _, isComplete) -> isComplete) states of
        (partial, [])           -> partial
        (partial, (complete:_)) -> partial ++ [complete]

Break the list of states, cutting before the first complete solution. Return the partial solution states that came before the cut, appending the complete solution state (the first after the cut if it exists). Ignore the rest.

ghci> printScan $ upToComplete $ scanAlgoXSimple' m1 []
([],
+---++---+---+---+---+---+---+---+
|   || 1 | 2 | 3 | 4 | 5 | 6 | 7 |
+===++===+===+===+===+===+===+===+
| 0 || 1 | 0 | 0 | 1 | 0 | 0 | 1 |
| 1 || 1 | 0 | 0 | 1 | 0 | 0 | 0 |
| 2 || 0 | 0 | 0 | 1 | 1 | 0 | 1 |
| 3 || 0 | 0 | 1 | 0 | 1 | 1 | 0 |
| 4 || 0 | 1 | 1 | 0 | 0 | 1 | 1 |
| 5 || 0 | 1 | 0 | 0 | 0 | 0 | 1 |
+---++---+---+---+---+---+---+---+
,False)

([0],
+---++---+---+---+---+
|   || 2 | 3 | 5 | 6 |
+===++===+===+===+===+
| 3 || 0 | 1 | 1 | 1 |
+---++---+---+---+---+
,False)

([3,0],
+--++---+
|  || 2 |
+==++===+
+--++---+
,False)

([1],
+---++---+---+---+---+---+
|   || 2 | 3 | 5 | 6 | 7 |
+===++===+===+===+===+===+
| 3 || 0 | 1 | 1 | 1 | 0 |
| 4 || 1 | 1 | 0 | 1 | 1 |
| 5 || 1 | 0 | 0 | 0 | 1 |
+---++---+---+---+---+---+
,False)

([3,1],
+---++---+---+
|   || 2 | 7 |
+===++===+===+
| 5 || 1 | 1 |
+---++---+---+
,False)

([5,3,1],
+--++--+
|  ||  |
+==++==+
+--++--+
,True)

This is more manageable. Recall the solutions are built backwards. So we tried row 0, then we tried row 3. That did not work out, so we backtracked and tried row 1, then row 3, then row 5, and that worked out.

Everything here is lazily evaluated, so scanAlgoXSimple' stopped searching when we stopped asking for further states after receiving the complete solution.

Here is a version of the function that only returns the solutions instead of all intermediate states.

algoXSimple' :: SparseMatrix -> [Key] -> [[Key]]
algoXSimple' (SparseMatrix rows activeCols) solution
    | size activeCols == 0 = [solution]
    | otherwise =
        [ s | (r, row) <- IntMap.toList rows,
              let m' = SparseMatrix (IntMap.filter (row `disjoint`) rows)
                                    (activeCols `difference` row),
              s <- algoXSimple' m' (r:solution) ]

ghci> printList "Set" sets1
+---++--------------------+
|   ||                Set |
+===++====================+
| 0 ||   fromList [1,4,7] |
| 1 ||     fromList [1,4] |
| 2 ||   fromList [4,5,7] |
| 3 ||   fromList [3,5,6] |
| 4 || fromList [2,3,6,7] |
| 5 ||     fromList [2,7] |
+---++--------------------+

ghci> algoXSimple' m1 []
[[5,3,1],[3,5,1],[5,1,3],[1,5,3],[3,1,5],[1,3,5]]

ghci> head $ algoXSimple' m1 []
[5,3,1]

Let's appreciate for a moment how far we have come so quickly and with what beautiful, terse code.

Compare this for compactness and simplicity with Dancing Links, which is how Algorithm X is usually implemented.

https://www.wikiwand.com/en/Dancing_Links

---

Adding column selection

Now, as Wikipedia says: to reduce the number of iterations, Knuth suggests that the column-choosing algorithm select a column with the smallest number of 1s in it.

I'm going to give this one all in one go:

(+.) = IntMap.unionWith (+)

algoX m = algoX' m []

algoX' :: SparseMatrix -> [Key] -> [[Key]]
algoX' (SparseMatrix rows activeCols) solution
    | size activeCols == 0 = [solution]
    | otherwise =
        let (c, colSum) = selectedColumn
        in  [ s | colSum > 0,
                  (r, row) <- IntMap.toList rows,
                  member c row,
                  let m' = SparseMatrix (IntMap.filter (row `disjoint`) rows)
                                        (activeCols `difference` row),
                  s <- algoX' m' (r:solution) ]
        where
            withOnes row'  = IntMap.fromSet (\_ -> 1) $ row' `intersection` activeCols
            colSums        = IntMap.foldl (\sums row'' -> sums +. withOnes row'')
                                          (IntMap.fromSet (\_ -> 0) activeCols)
                                          rows
            selectedColumn :: (Key, Int)
            selectedColumn = snd . minimum $ map (\(key, cSum) -> (cSum, (key, cSum)))
                                                 (IntMap.toList colSums)

Playing in the interpreter with some of the new expressions is left as an exercise to the reader.

This version is a little verbose with selecting the column, but the main part of the algorithm remains clean. We've added guards: if the column sum is zero or if the selected column is not in the row then we cut that search path. (We could be caching some things to make the column selection computation more efficient...)

ghci> head $ algoX m1
[5,3,1]

Taking the first solution, it is what we expected.

One thing to note is that because we are now selecting one column for any particular matrix, we do not search through all permutations of our single solution anymore:

ghci> algoX m1
[[5,3,1]]

If there is interest, we can throw these algorithms at some large Sudoku-style puzzles in some future post.

---

Play around yourself

The code is available here:

https://github.com/AlexeyFeigin/algorithm-x

The Demo module includes all the variables defined in this blog post.

You will need to install ghc (the Glasgow Haskell Compiler) and cabal (the Haskell package manager), which both come with the Haskell Platform:

https://www.haskell.org/platform/

Happy hacking!

---

_{Attribution: Cover photo by CHUTTERSNAP.}

---

Appendix

Printing functions

(!.) :: Num a => Row -> Key -> a
row !. key
    | member key row = 1
    | otherwise      = 0

toFilledList :: ActiveCols -> Row -> [Int]
toFilledList activeCols row = map (row !.) $ elems activeCols

instance Show SparseMatrix where
    show (SparseMatrix rows activeCols) =
        ('\n' :) $ render show show show $ Table (Group NoLine     $ map Header $ IntMap.keys rows)
                                                 (Group SingleLine $ map Header $ elems activeCols)
                                                 (map (toFilledList activeCols) $ IntMap.elems rows)

printRow :: Row -> IO ()
printRow row = print $ SparseMatrix (IntMap.fromList [(0, row)]) row

oneColumnTable :: String -> [(key, value)] -> Table key String value
oneColumnTable columnHeader ls = Table (Group NoLine $ map (Header . fst) ls)
                                       (Header columnHeader)
                                       (map (\(_, s) -> [s]) ls)

printList :: Show a => String -> [a] -> IO ()
printList title = putStrLn . render show id show . oneColumnTable title . zip [0::Int ..]

printScan :: Show a => [a] -> IO ()
printScan = putStrLn . intercalate "\n\n" . map show

Updates

• 2021-12-17 Tweaked tags