After inventing a new shinier sort of wheel and writing my own string class in C++, it was time to re-implement some OR algorithms, and I spent some time creating F# code for the Transportation Problem. The calculations are mostly routine and not very interesting, but at the heart of it all is a little mechanism that I found lends itself to be coded in a very succinct and elegant way, with a few standard library calls and little more; I owe to Microsoft and the F# team my very knowledge of the word succinct, and it is one of the three main characteristics of the language (the others being expressive and efficient, but I was familiar with those).
I may or may not show the complete code for the transportation algorithm on a later occasion, but for now I am just pulling this little mechanism out of its context and abstract it for isolated treatment: the Stepping Stone algorithm.
Recap
Very quickly, because it is not necessary in order to understand the code: In the transportation problem we have a number of suppliers S1,..,Sn that each supply a certain quantity of something. And we have a number of destinations or sinks D1,…,Dm that consume certain amounts of said something. m and n can and will be different. The something needs to be transported from the Si to the Dj (i in 1…m, j in 1…n) and each Si can deliver to more than one Dj. Equally, each Dj can receive from more than one Si. The total amount of something delivered is the same as the total amount received. Delivering from a certain Si to a Dj costs Cij per unit of something.
The problem: Find out how much each Si has to deliver to each Dj so that everyone is happy and the total cost is minimal.
The way this is done:
- Find some solution that satisfies everyone, without necessarily having minimal cost. This is easy.
- Check if this is optimal. This needs a bit of calculation.
- If not optimal, find a new ‘satisfactory’ solution which is a bit better, by switching one supplier/destination Si-Dj against another Sk-Dl. Then rinse and repeat.
Step 3 is the one that is of interest here.
Representation
The problem can be solved as a normal linear program using the simplex algorithm. However, the coefficient matrix for this has a very special form that allows for a more compact (or even succinct) representation, as a transportation tableau.
Fig. 1
Each row in Fig. 1 represents all the deliveries of one supplier, and each column is what one particular sink consumes. Normally, this would be where we see the cost for each combination of the two, but this is not really needed at this point.
The orange cells represent a first solution to the problem, i.e. combinations of suppliers and destinations where everyone is satisfied. However, at this point the transportation algorithm may want to switch to a different solution by, say, having supplier S3 deliver to destination D1 (and you will have noticed how here in programming world we start counting from 0 and not from 1). This basis change happens by finding a path from the new (entering) cell through all the cells in the current solution and back to the new cell, with a few conditions. Fig. 2 shows the tableau with the new cell in blue.
Fig. 2
Here is how that path needs to be constructed:
1. Start from one entering cell and find another cell in the same row or column that is in the solution (orange).
2. Continuing from the selected cell, find the next one within the solution (orange). If you moved within a row to find the current cell, now you need to check within that cell’s column; if you found it in the same column as the one before, now the next one must be found in the same row as the current one.
If in the process you reach the entering cell again, then you’re done. Otherwise keep repeating step 2. The path we needed to find for entering cell (3,1) is shown in Fig. 3.
Fig. 3
The Mathematics here is not quite complete – we are not (yet) solving the transportation problem as such, don’t make use of the cost matrix, or of the supplier and sink capacities, and are only concentrating on a particular step of the algorithm. But this can be nicely isolated and formulated as a coding problem without all the additional information. Once this path finding operation is out of the way, the remaining calculations are rather trivial, and we might fill them in later.
The Coding Task
Given a list of solution cell co-ordinates of the form (row, column), and a single co-ordinate for the entering cell that can not be part of the solution, find a valid path from the entering cell back to itself through any number of solution cells, with the following restrictions:
- Only move horizontally or vertically (no diagonal movements), and
- After moving along a row, you have to move along a column, and vice versa.
Our input data looks as follows:
let solution = [(0,0);(0,1);(2,0);(2,2);(2,3);(1,2);(3,3)]
let entering = (3, 1)
Data Structures
While performing the algorithm, the code needs to search the tableau by row and by column, potentially many times, so the representation as a list is not ideal. It would require a linear search each time we need to pull out a row or a column. This could be made a bit better by sorting it first, but that could only work for either rows or columns. If, on the other hand, we create one big matrix or 2-dimensional array, the row and column splicing would be pretty easy and efficient – however, this would waste lots of space. In a real transportation problem, we have
n = number of suppliers = number of rows
m = number of destinations = number of columns
m+n-1 = number of cells in any solution
m*n = number of cells in the matrix
(m+n-1)/(m*n) = occupation of the matrix
The bigger the problem gets, the less space of the total matrix is actually needed.
4 suppliers, 4 destinations –> 43.75%
10 suppliers, 20 destinations –> 14.5%
100 suppliers, 200 destinations –> 1.495%
This is a case where at first glance a sparse matrix implementation would come in handy. However, they have their own issues, most importantly in this case the lack of efficient row and column splicing (they tend to have one of the two at most). We create our own problem-specific implementation:
- One lookup table (dictionary) with a list of column indices for each row index
- Another lookup table (dictionary) with a list of row indices for each column index
Creating the Dictionaries
What we have is
let solution = [(0,0);(0,1);(2,0);(2,2);(2,3);(1,2);(3,3)]
and, just considering one of the two cases, what we want is something that looks similar to
(* pseudo code *)
let (lookuptable: keyvaluetype) = [(0,[0;1]);(1,[2]);(2,[0;2;3]);(3,[3])]
Luckily, this can be done by piping the input list through a few standard library functions. First, let F# group the tuples of the input list by some key. The Seq.groupBy functions wants a function argument that calculates for each element what its key value is, and then does the grouping. Just looking at our first case (“One lookup table (dictionary) with a list of column indices for each row index”) the key to group by has to be the row index, so the function we provide to Seq.groupBy pulls the row index out of each tuple.
solution |> Seq.groupBy (fun (r,_) -> r)
And, remembering that such a function already exists in the library, this can be abbreviated to
solution |> Seq.groupBy fst
This creates a collection of all elements that correspond to each row index. Since the elements are the (row, column) tuples, what we have now is
(* pseudo code *)
[(0,[(0,0);(0,1)]);(1,[(1,2)]);(2,[(2,0);(2,2);(2,3)]);(3,[(3,3)])]
It is close to, but not quite, what we want. Seq.map is used to convert a sequence of some sort of elements into a sequence of some other sort, so this is what we need. Make the sequence of tuples where the second element is a list of index tuples into one where the second element is a list of indices.
(* ... *)
|> Seq.map (fun (key, valSeq) ->
(key, valSeq |> Seq.map (fun(_,c) -> c)))
We are using Seq.map twice, once to transform the ‘outer’ sequence, and then for the sequence that is contained in each of the tuples that are the elements. Again, there’s a library replacement for the function that pulls out the column indices and the code can be shortened to:
(* ... *)
|> Seq.map (fun (key, valSeq) ->
(key, valSeq |> Seq.map snd))
The new collection is
(* pseudo code *)
[(0,[0;1]);(1,[2]);(2,[0;2;3]);(3,[3])]
In the end we can use the dict function to create a lookup table out of this. A complete function that takes a solution list as in our input data and transforms it into a lookup table like those that we decided we need for the algorithm, could be implemented like this:
let makeDict lst =
lst
|> Seq.groupBy fst
|> Seq.map (fun (key, valSeq) -> (key, valSeq |> Seq.map snd))
|> dict
Now we need to do the same, but for the indices being on the columns. The function would look exactly the same, only for the roles of the calls to fst and snd swapped. A more generic version of the function is in order, where this information can be passed into instead of being hard coded:
(* The more general function *)
let makeDict keySelector valSelector lst =
lst
|> Seq.groupBy keySelector
|> Seq.map (fun (key, valSeq) -> (key, valSeq |> Seq.map valSelector))
|> dict
Path Finding
Now that we can look up any column for a row index and any row for a column index, it is pretty straightforward to take our entering cell and search for a path back.
let findPath (targetRow, targetCol) lstSolution =
(* Create the lookup tables for rows and columns.
Both use makeDict but pass in the 'selector'
functions in different order *)
let rowDict = lstSolution |> makeDict fst snd
let colDict = lstSolution |> makeDict snd fst
(* From current position, find the row we're in,
then go through all cells (except our current one).
If we end up on the same column as the entering
cell, aka 'target', then this is the last element
of the list we are creating. Otherwise, check to
see if we can move on in this column. If both fails,
there's nothing there there. *)
let rec findInRow currentRow currentCol =
rowDict.[currentRow]
|> Seq.tryPick (fun col ->
if col=currentCol then None
else if col=targetCol then Some([(currentRow, col)])
else
match findInCol currentRow col with
| Some(l) -> Some((currentRow, col)::l)
| _ -> None)
(* Like findInRow, but looking through current column *)
and findInCol currentRow currentCol =
colDict.[currentCol]
|> Seq.tryPick (fun row ->
if row=currentRow then None
else if row=targetRow then Some([(row, currentCol)])
else
match findInRow row currentCol with
| Some(l) -> Some((row, currentCol)::l)
| _ -> None)
(* Return value is created by calling one of the
nested functions. *)
findInCol targetRow targetCol
This uses the Seq.tryPick function. Seq.tryPick iterates through a sequence and calls a specified function for each element. That function calculates a value from its argument and wraps it up in an option type value, or returns None. For None, Seq.tryPick then goes to the next element of the sequence, otherwise the returned option type value becomes the return value of Seq.tryPick. This is very similar to Seq.tryFind, except that Seq.tryFind uses a predicate function just to check if a sequence element is the one to be picked, and then wraps this element itself into an option type value; Seq.tryPick is a bit more flexible in that it allows to return something that is not the element value itself.
The thing we want to return here is a list with the elements of the found path. When we successfully arrive at the end of a valid path, the current cell becomes the end of that list, in
else if col=targetCol then Some([(currentRow, col)])
or
else if row=targetRow then Some([(row, currentCol)])
If the current cell is not the end of the path, but if a valid path is found ‘from here’, then it is prepended to the list that has been returned in the recursive foundInRow or foundInCol call.
match findInRow row currentCol with
| Some(l) -> Some((row, currentCol)::l)
| _ -> None)
If this recursive call has not found anything, then the current cell is not part of the path, and returns None to inform its own caller accordingly. The top-level caller of this will either return an option type value with a list of cells that show the valid path, or None. It may be instructive to see how the algorithm finds its path by using backtracking: the way the recursive calls are made results in a trial-and-error movement along the cells in the current solution.
Backtracking
Fig. 4 shows how, starting from entering cell (3,1) in step [1], the algorithm looks through the current column and finds a cell in step [2], then looks in that row [3], then in the new column [4], then in the row and tries the first cell it finds there [5] and looks through that cell’s column [6]. The path traced so far is represented in the state of the code through the current call stack of the recursive calls to findInRow and findInCol. However, from the last cell that was found (1,2), there is no way forward.
Step [6] did not lead to a result (see Fig. 5) and another try has to be made starting from [5]. However, there’s nothing else to try in that column, so step [5] is itself invalidated. The process has to re-start from [4], and as Fig 6 shows, new steps [5] and [6] are found that result in the required path back to the entering cell.
It can be instructive to see the backtracking process unfolding during program execution. The following is a modified version of the findPath function, that is functionally identical, but prints out its current activity.
let findPath (targetRow, targetCol) lstStones =
let rowDict = lstStones |> makeDict fst snd
let colDict = lstStones |> makeDict snd fst
let rec findInRow tab currentRow currentCol =
rowDict.[currentRow]
|> Seq.tryPick (fun col ->
if col=currentCol then None
else
printfn "%s%A" tab (currentRow, col)
if col=targetCol then
printfn "%sFound." tab
Some([(currentRow, col)])
else
match findInCol (tab + "..") currentRow col with
| Some(l) -> Some((currentRow, col)::l)
| _ -> None)
and findInCol tab currentRow currentCol =
colDict.[currentCol]
|> Seq.tryPick (fun row ->
if row=currentRow then None
else
printfn "%s%A" tab (row, currentCol)
if row=targetRow then
printfn "%sFound." tab
Some([(row, currentCol)])
else
match findInRow (tab + "..") row currentCol with
| Some(l) -> Some((row, currentCol)::l)
| _ -> None)
findInCol "" targetRow targetCol
Example
Fig. 7 
Fig. 8
The current solution as shown in Fig. 7 is quite realistic, for a small problem size, and it would look like this after 2 or 3 steps into the transportation algorithm. While it seems quite benign, we can see from Fig. 8 that the path belonging to an entering cell can become a bit complex. Here is the data and code to try out this example (assuming that the earlier mentioned functions are implemented somewhere.
let current = [(1, 3); (1, 0); (2, 0); (2, 1); (10, 1); (10, 5); (9, 5);
(9, 4); (7, 4); (5, 3); (5, 2); (3, 2); (0, 0); (4, 2);
(6, 3); (8, 4); (11, 5); (11, 6); (12, 6); (13, 6);
(13, 7); (14, 7); (14, 8); (15, 8); (16, 8); (16, 9);
(17, 9); (18, 9); (19, 9)]
let entering = (15, 3)
current |> findPath entering |> printfn "%A"
Remarks
- This is only a part of the overall transportation algorithm, which I think merits some exploration on its own. I will post the whole thing when I have some more time. (Update 18-Jan-2011: This has now been done, see Transportation Algorithm in F#)
- Instead of using the self-made lookup tables, this really shouts for the use of some standard sparse matrix implementation. This will be topic of a future post.
- In the real transportation problem, if it is balanced and solutions are not degenerate, there exists one unique path for the entering cell, meaning that a path can always be found and that it doesn’t matter if you start looking in the row or column. In arbitrary examples this may be different, so you can have situations where there is no path, or where the one you find is only one possibility, and where there may be a difference between kicking this off with a call to findInRow or to findInCol.
- There are a couple of alternative implementations for the findPath function, which I talk about here.
Fig. 4
Fig. 5
Fig. 6
