Dynamic Programming

I've started reading Jones & Pevzner, An Introduction to Bioinformatics Algorithms.

These authors spend substantial time on a classic computer science method called "dynamic programming" (invented by Richard Bellman). It is widely used in bioinformatics. For example, sequence alignment algorithms such as Needleman-Wunsch and Smith-Waterman are dynamic programming methods. The key is to have a problem that can be broken into stages or steps; the problem is attacked from one end or the other and the result for each stage is computed in turn, based on current knowledge, and the result is stored.

A problem they use to illustrate the method is the "change" problem, which I first ran across in SICP. (Another is Towers of Hanoi).

As much as I enjoyed the example, I think it is an unfortunate choice for the biologically oriented student of bioinformatics, because it is really more complex than necessary to illustrate the dynamic programming idea; the complexities are not relevant to the sequence problem, and may well distract the naive student.

I would talk first about a simpler problem. Imagine a triangle made up of elements with integer values, having one element on the first line, two on the second, three on the third and so on. Given a current position, we must move down and either left or right to the next element.


Our problem statement is to evaluate alternative paths from top to bottom through this triangle (e.g. to discover the path with the maximum total value), as described in this problem from Project Euler.

As you can appreciate, triangles can become so large that it is impossible to enumerate and evaluate all paths. The number of paths is 2**(n-1) for a triangle with n levels, and this value gets large rapidly with increasing n. To solve the problem, we work from the bottom to the top.

However, what I really want to talk about in this post is the change problem, which is a version of something more generally known as the knapsack problem.

We have a container (the knapsack) with a known capacity, and objects of various size which we wish to store in it. We may desire to know different things, for example, how many different combinations of objects can we find that completely fill the knapsack? This is our problem statement for the change problem: given American coins (1, 5, 10, and 25 cents), how many ways are there to make change for a dollar?

Rather than go back and look up the answer in SICP (it's been a long time, and I was never much good with Scheme), I decided it would be more fun to try to solve it without looking at any other resources. The key for me was to visualize the problem as a tree, in which the top-level node contains a set of 100 pennies. What can we do? Well, we could change 5 pennies to a nickel, or 10 pennies to a dime, and so on. Each one of these new forms of change for a dollar becomes a sub-node below the first one. From each sub-node, in turn, we attempt further transformations, which will become sub-nodes in turn.

For example, starting with only 10 pennies (P=penny,N=nickel,D=dime):

10P 0N 0D => 5P 1N 0D => 0P 2N 0D => 0P 0N 1D => 0P ON 1D

Two things to notice: we can do some kind of depth first traversal, or bread-first. (My graph algorithm skills are weak, but that much seems clear at least). And the depth-first would presumably use recursion (for which this problem is famous). The other thing to notice is that we can get duplicates---there are two paths to the 1 dime state.

My solution (Python code) follows below. Simple things: I use a dict to hold each single kind of change, and a utility function that checks whether it is possible for us to change, say, pennies to a nickel for a current dict. This function returns either a new dict containing the result or None. The loop which solves the problem starts with a list containing a single element, a dict like this:

{ 1:100,5:0,10:0,25:0 }

The code does the following:

for each dict that we have discovered but not tried to manipulate: save it in our permanent list of solutions work with it further... for each coin type i from penny to dime: for each coin type j from the next highest to the last one: try to make change from i to j if we succeed: check to be sure this is not a solution we know already add it to the list for the next round and mark "not done"

Notably, this solution adaptable to other coin types and other total amounts. However, it seems not to be very efficient, since it slows dramatically with amounts beyond a few dollars. I'll have to work on that.

# how many ways to change a dollar? # each possible combination is a dict coins = [1,5,10,25] N = 100 L = [ { 1:N,5:0,10:0,25:0 } ] # make substitution if possible # change n * coin 1 => coin 2 def change(D,c1,c2): rD = dict() for k in D: rD[k] = D[k] # 10 => 25 is special if c1 is 10 and c2 is 25: if rD[10] < 2 or rD[5] < 1: return None else: rD[5] -= 1; rD[10] -= 2 rD[25] += 1 else: if not c1*rD[c1] >= c2: return None rD[c1] = rD[c1] - c2/c1 rD[c2] += 1 return rD # - - - - - - - - - - - - - - - - - - - - - - - - - - # start with N pennies # for each coin 1..10, try changing to higher denom # save results in L to explore further changes... # keep a temp copy separate from L so we can append in loop temp = list() results = list() while True: #print 'working:', str(len(L)).rjust(5) done = True for D in L: # we already checked, D is not yet in results results.append(D) # for each possible change given coins available for i in range(len(coins)-1): for j in range(i,len(coins)): # change if possible rD = change(D,coins[i],coins[j]) #if we haven't seen this one yet if rD and not (rD in results or rD in temp): temp.append(rD) # since there is a new one(s) do not quit done = False if done: break L = temp[:] temp = list() print 'finished', str(len(results)).rjust(5) results.sort(reverse=True) for D in results: print D

finished 242 {1: 100, 10: 0, 5: 0, 25: 0} {1: 95, 10: 0, 5: 1, 25: 0} {1: 90, 10: 0, 5: 2, 25: 0} ... {1: 0, 10: 5, 5: 0, 25: 2} {1: 0, 10: 0, 5: 0, 25: 4}