Learning F#: Case study with branch and bound – Part I (page 3)

But before we get back to coding, let’s understand a bit more about branch and bound. In the figure, we have a three-variable problem (binary in that each variable has two set positions, for our 0-1 Knapsack these are 0 and 1 of course). Application of branch and bound to this problem dates back many years. In Part II, we will develop branch and bound more generally but applying it first to the 0-1 knapsack will aid our incremental development in F#.

Branch and Bound Illustration

Suppose that the function g can tell you (at least estimate) what the best possible solution (max in our case) might be for any solution combination of variables given the ones that are set already. For example, in the figure we see that the value of g once we have selected v0 to be 0 is G=79. Suppose further that via our exploration of the variable settings we have determined a solution with utility of U=90 (following the red line here, this would be for settings of v0 = 1, v1 = 0, and v2 = 0). Since we have already found a solution with a utility of 90, there is no reason to explore the branch starting at the node where g is 79 as there are no solutions there that can be better than 79. This assumes of course that our g estimate is good, namely that it doesn’t underestimate the potential solutions for any subsequent variable settings. The tighter g is to the actual best solution within its sub-space the better, but the method works as long as the estimator does not underestimate.

This is the basis for search savings in branch and bound, if we have estimators that can give what the best possible solution is for subsets of variables (at least an estimator that doesn’t underestimate possible solutions) then once we start finding actual solutions we don’t need to explore spaces that cannot be any better. As reference, this web video walks through the branch and bound procedure for 0-1 knapsack; we will move on with coding.

This wiki page gives a generic (pseudo code) version of branch and bound, although this version assumes that we want to minimize the objective function not maximize as we are for our knapsack. We will see later that handling either maximizing or minimizing is not a problem, but for now let’s focus on maximizing.

Utility estimate of a knapsack with partial variable settings

So in order for our branch and bound routine to work, we need to estimate the best value that a partial variable setting would have. In the figure above, this would be the ‘G’ values at each node. Recall that each node is a setting of variables such that each has one of the Unset, Zero or One value. In the example, we somehow determined that after setting v0 = 0 that G=79 (this is an arbitrary number for demonstration). How can we get such an estimate for our problem? We already have a function multiplySumUtility that computes the utility of a knapsack given selections, and a function that determines the weight of such as well.

Without getting into the details, our estimator can work like this:

  1. Sort the variables in order of best utility to weight ratio
  2. Determine the utility of the variables that have already been assigned, call this assignedUtility
  3. Determine the weight of the assigned variables, and set remainingWeight to the weight limit less this weight
  4. For the Unset variables var (in order according to step 1), and while the remainingWeight is > 0:
    • if var.Weight <= remainingWeight, add var.Weight to unassignedUtilityEstimate and reduce remainingWeight by this weight
    • otherwise take a fraction of this items value: let weightFrac = remainingWeight / var.Weight; add weightFrac*var.Utility to unassignedUtilityEstimate; set remainingWeight = 0

There are other forms of knapsack, like fractional and unbounded, but aren’t considering them now

This estimates the utility of a node as: the sum of the existing utility for variables that have been assigned plus the best possible value of those that have not. This is done by adding the best possible choices in terms of utility to weight ratio including a fractional one at the end. Of course in our real knapsack we can only choose or not choose an item, but this works as our estimator, which cannot be less that the actual utility for the node. That is, for these settings this is the best that could possibly be obtained for any sub-node in the search tree.

So let’s consider the F# coding for this:

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 

let estimator wlimit vars =
    // assume vars are sorted in decreasing utility/weight ratio
    let assignedUtility = multiplySumUtility vars
    let mutable unassignedUtilityEstimate = 0.0
    let mutable remainingWeight = wlimit - multiplySumWeight vars
    Array.iter (fun var -> if remainingWeight > 0.0 then do
                             if var.Weight <= remainingWeight then do
                               remainingWeight <- remainingWeight - var.Weight
                               unassignedUtilityEstimate <- unassignedUtilityEstimate + var.Utility
                             else do
                               let frac = var.Weight/remainingWeight
                               unassignedUtilityEstimate <- unassignedUtilityEstimate + frac * var.Utility
                               remainingWeight <- 0.0
                             ) vars

    assignedUtility + unassignedUtilityEstimate

This seems to be on track, we basically use the Array.iter to go through each var, in decreasing order of utility to weight ratio (assumed), an adjust the mutable values. But, this doesn’t work — basically saying that mutable variables cannot be used in a closure. Well, the closure is the unnamed function we used as the Array.iter argument and of course the mutable variables are those we defined. We can fix this using ref variables instead as so:

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 

let estimator wlimit vars =
    // assume vars are sorted in decreasing utility/weight ratio
    let assignedUtility = multiplySumUtility vars // value from set vars
    let unassignedUtilityEstimate = ref 0.0
    let remainingWeight = ref (wlimit - multiplySumWeight vars)
    Array.iter (fun var -> if !remainingWeight > 0.0 && var.Setting = Unset then do
                             if var.Weight <= !remainingWeight then do
                               remainingWeight := !remainingWeight - var.Weight
                               unassignedUtilityEstimate := !unassignedUtilityEstimate + var.Utility
                             else do
                               let frac = var.Weight / !remainingWeight
                               unassignedUtilityEstimate := !unassignedUtilityEstimate + frac * var.Utility
                               remainingWeight := 0.0
                             ) vars

    assignedUtility + !unassignedUtilityEstimate

Ref variables can be used in closures, recall these are assigned using the := and accessed using the ! operator as seen. But is this the right way to go? Idiomatic F# tries to stay away from ref variables (although of course sometimes they are required and the best way to do things), so if we need ref variables we should at least think about other ways to do what we need. What if we used an Array.fold to ‘accumulate’ a result across the elements, like:

Is this version clearer in terms of easily seeing what the routine does?

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 

let estimator wlimit vars =
    // assume vars are sorted in decreasing utility/weight ratio
    let assignedUtility = multiplySumUtility vars
    let remainingWeight = wlimit - multiplySumWeight vars
    let (leftOverWeight,unassignedUtilityEstimate) = 
        Array.fold (fun acc var -> match acc,var.Setting with
                                   | (_,accUtil),Zero   // we can't use _,accUtil,Zero here - that's a 3-tuple 
                                   | (_,accUtil),One
                                   | (0.0,accUtil),_ -> (fst acc, accUtil)
                                   | (accWt,accUtil),varSetting -> 
                                       if accWt >= var.Weight then
                                          (accWt - var.Weight, accUtil + var.Utility)
                                       else 
                                          let frac = var.Weight / accWt
                                          (0.0, accUtil + frac * var.Utility)
                   ) (remainingWeight,0.0) vars

    assignedUtility + unassignedUtilityEstimate

F# note: The fst function returns the first element of a tuple

The Array.fold function uses a base element, here the tuple (remainingWeight,0.0) and applies the given function to each element threading an accumulator through the result (called acc here) [yes, this gives the same result as before]. The tuple is a pair such that the first part is the remainingWeight (or ‘leftOverWeight’ as it adjusts) and the second is the estimated value. Now our function uses pattern matching to determine the next accumulator return value. The idea of the pattern is: if the variable is not Unset (e.g. Zero or One) or if the accumulator’s first value (the leftOverWeight) is zero, then just pass along the acc as is, we don’t want to change it — we are only adding contributions here from un-set variables. Otherwise, apply our calculation for adjusting each of the leftOverWeight and estimated utility value.

F# note: the order of pattern matching matters

Notice that in the pattern match we had to be careful to match each of the acc and var.Setting elements properly. For example, the parens in (_,accUtil),Zero are needed so that the parenthetical part matches acc and the Zero matches var.Setting, if we tried to use _,accUtil,Zero it would be an error since that is a 3-tuple and not a tuple of a 2-tuple and a DiscreteVar variable — which is what acc,var.Setting is.

0-1 Knapsack branch and bound

Finally, with our estimator in place and some adjustments to the exhuastive search code we’ll do here, we can tackle branch and bound.

F# note: Here is a handy operator reference

We still need to determine variables that are unset, so for the purposes of covering more F#, let’s revisit the unsetVar function from before. That returned a tuple of ‘any-unset, index-of-unset’ basically. As was hinted, instead let’s just return an int such that if it is -1 then we interpret that as no vars are unset. And let’s explore using Array.findIndex, but make sure we capture the exception thrown:

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 

open System.Collections.Generic

(* OTHER CODE *)

(* return index of unset (first) var, or -1 if all are set (to One or Zero) *)
let unsetVar vars =
    try
        Array.findIndex (fun var -> var.Setting = Unset) vars
    with
        | :? KeyNotFoundException -> -1

As in other computer languages, F# uses a try/catch approach (using keywords try and with however). The new unsetVar routine uses the Array.findIndex with a predicate (condition) of looking for an Unset variable. It will throw an exception of KeyNotFoundException if not found which is caught in the with clause. Note that the System.Collections.Generic library must be included (this is done with the open keyword) in order to bring the KeyNotFoundException name into scope. The :? operator returns true if the argument matches and otherwise false, so if a KeyNotFoundException exception is raised it will be caught. Our new unsetVar routine is simpler and problably more readily understood than before.

Since we used a queueing style, we can readily adapt our exhaustive searcher to a branch and bound solver. We just need to incorporate the estimator function such that it prevents partial solutions that are already worse than a solution we have from going on the queue. Here it is, with our new unsetVar function.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
48: 
49: 
50: 
51: 
52: 
53: 
54: 
55: 
56: 
57: 
58: 
59: 
60: 
61: 

let knapsack01 wlimit vars initLowestValue =
// make sure that vars are sorted in utility/weight descending order [Note: minus sign]
let varsSorted = Array.sortBy (fun elem -> -elem.Utility / elem.Weight) vars   
// same as exhaustive
let mutable bestUtility = initLowestValue;
let mutable queue = [varsSorted]
let mutable bestVars = [| |]     // best var settings found, nothing so far.
// these are only Informational, not needed in the routine
let mutable numberOfQueues = 1   // we just queued the start, so 1
let mutable numberOfEvals = 0    // number of times we evaluated the function with all vars set
(* return index of unset (first) var, or -1 if all are set (to One or Zero) *)
let unsetVar vars =
try
Array.findIndex (fun var -> var.Setting = Unset) vars
with
| :? KeyNotFoundException -> -1
while not queue.IsEmpty do
match queue with
| curVars :: restOfQueue -> 
let varIndex = unsetVar curVars 
if not (varIndex = -1) then
let mutable newQueue = restOfQueue  // if the 'g' funcs pass later, we'll add to the queue, 
// otherwise its just the rest of the queue
let tvar = Array.get curVars varIndex
// change the unset var to One and Zero
let zeroVar = { tvar with Setting = Zero }
let varArrayZero = Array.copy curVars  // copy the array so we can change element
Array.set varArrayZero varIndex zeroVar
let gEstimateZero = estimator wlimit varArrayZero
if gEstimateZero > bestUtility then 
newQueue <- varArrayZero :: newQueue
numberOfQueues <- numberOfQueues + 1
let oneVar = { tvar with Setting = One }
let varArrayOne = Array.copy curVars  // copy the array so we can change element
Array.set varArrayOne varIndex oneVar
let gEstimateOne = estimator wlimit varArrayOne
if gEstimateOne > bestUtility then 
newQueue <- varArrayOne :: newQueue
numberOfQueues <- numberOfQueues + 1
queue <- newQueue   // set the queue to the new queue
else
let util = multiplySumUtility curVars
let wt = multiplySumWeight curVars
numberOfEvals <- numberOfEvals + 1
if util > bestUtility && wt <= wlimit then
bestUtility <- util
bestVars <- curVars
queue <- restOfQueue // we've processed curVars, take it off the queue
| _ -> ()  // queue isempty, do nothing
printfn "Total settings queued = %d, total number of evals: %d" numberOfQueues numberOfEvals
(bestUtility, bestVars)

Running this, we get the same solution as our exhaustive solver. However, now we’ve only visited 16 nodes (versus 31 before) and needed to evaluate 4 full settings (versus all 16 combinations before).

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 

Total settings queued = 16, total number of evals: 4
Best Utility = 14.000000
{Name = "map";
Weight = 0.5;
Utility = 3.0;
Setting = One;}
{Name = "gps";
Weight = 1.0;
Utility = 3.0;
Setting = One;}
{Name = "food";
Weight = 5.0;
Utility = 8.0;
Setting = One;}
{Name = "tent";
Weight = 14.5;
Utility = 5.0;
Setting = Zero;}

Our code is getting a bit complex and probably should be tested and timed for improvements, but that is a subject for Part III. Let’s take a look at other improvements. The code concerning adapting the queue is a bit convoluted, are there better ways to handle this? We have seen pattern matching against discriminated unions (to determine which setting applies), and for given values, but there is another option active patterns. An active pattern uses code execution to determine which case applies.

Take a look at the code on lines 22-30 below. This is an active pattern definition that takes a 3-tuple as an argument and matches to four cases, TakeBoth, TakeZero, TakeOne or TakeNeither. The if expression decides which pattern to match based on the argument. In this case, for example, if both gZero (this will be the estimate for the Zero case in the body of the code) and gOne are greater than the bestSoFar then TakeBoth is matched. A similar explanation applies to the other cases, and the discussion of its application follows.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
48: 
49: 
50: 
51: 
52: 
53: 
54: 
55: 
56: 
57: 
58: 
59: 
60: 
61: 
62: 
63: 
64: 
65: 
66: 
67: 
68: 
69: 
70: 
71: 
72: 
73: 
74: 
75: 
76: 

let knapsack01 wlimit vars initLowestValue =
// make sure that vars are sorted in utility/weight descending order [Note: minus sign]
let varsSorted = Array.sortBy (fun elem -> -elem.Utility / elem.Weight) vars   
// same as exhaustive
let mutable bestUtility = initLowestValue;
let mutable queue = [varsSorted]
let mutable bestVars = [| |]     // best var settings found, nothing so far.
// these are only Informational, not needed in the routine
let mutable numberOfQueues = 1   // we just queued the start, so 1
let mutable numberOfEvals = 0    // number of times we evaluated the function with all vars set
(* return index of unset (first) var, or -1 if all are set (to One or Zero) *)
let unsetVar vars =
try
Array.findIndex (fun var -> var.Setting = Unset) vars
with
| :? KeyNotFoundException -> -1
(* active patterns *)
let (|TakeBoth|TakeZeroOnly|TakeOneOnly|TakeNeither|) (gZero, gOne, bestSoFar) = 
if gZero > bestSoFar && gOne > bestSoFar then
TakeBoth
else if gZero > bestSoFar then
TakeZeroOnly
else if gOne > bestSoFar then
TakeOneOnly
else
TakeNeither
while not queue.IsEmpty do
match queue with
| curVars :: restOfQueue -> 
let varIndex = unsetVar curVars 
if not (varIndex = -1) then
let tvar = Array.get curVars varIndex
// change the unset var to One and Zero
let zeroVar = { tvar with Setting = Zero }
let varArrayZero = Array.copy curVars  // copy the array so we can change element
Array.set varArrayZero varIndex zeroVar
let gEstimateZero = estimator wlimit varArrayZero
let oneVar = { tvar with Setting = One }
let varArrayOne = Array.copy curVars  // copy the array so we can change element
Array.set varArrayOne varIndex oneVar
let gEstimateOne = estimator wlimit varArrayOne
// update the queue
queue <- match (gEstimateZero, gEstimateOne, bestUtility) with
| TakeBoth -> 
numberOfQueues <- numberOfQueues + 2
varArrayOne :: varArrayZero :: restOfQueue
| TakeZeroOnly -> 
numberOfQueues <- numberOfQueues + 1
varArrayZero :: restOfQueue
| TakeOneOnly -> 
numberOfQueues <- numberOfQueues + 1
varArrayOne :: restOfQueue
| _ -> restOfQueue // Neither case
else
let util = multiplySumUtility curVars
let wt = multiplySumWeight curVars
numberOfEvals <- numberOfEvals + 1
if util > bestUtility && wt <= wlimit then
bestUtility <- util
bestVars <- curVars
queue <- restOfQueue // we've processed curVars, take it off the queue
| _ -> ()  // queue isempty, do nothing
printfn "Total settings queued = %d, total number of evals: %d" numberOfQueues numberOfEvals
(bestUtility, bestVars)

??? I think the active pattern version might be a bit clearer, but it relies on also understanding the active pattern. Is there still a better way?

Lines 51-61 show use if the active pattern, here the queue is updated depending on which case (active pattern) applies. Perhaps this is more readily understandable than the previous approach?

Also to circle back to an earlier decision, we made the Setting field of DiscreteVar mutable with the idea that we would need to change it, but we never did. Note that we create copies of an unset var (called tvar in the code) and use the with clause to change the setting to a new DiscreteVar. So we never actually update the Setting in place, so here is the minor update:

1: 
2: 
3: 
4: 
5: 
6: 
7: 

type DiscreteVar = 
{                               (* this is a record type in F# *)
Name : string;
Weight : float;
Utility : float;
Setting : ZeroOneVarSetting;         // discrete value setting
}