January 2018

S M T W T F S
  123456
78910111213
14151617181920
21222324252627
28293031   

Style Credit

Expand Cut Tags

No cut tags
Thursday, August 26th, 2010 03:32 pm
A couple of days ago, Chris Yocum asked me to help him out with some OCaml code. He's learning OCaml at the moment, and thought he'd practice it by solving a problem he'd encountered in his research on early Irish law.

In early Irish law, murder was punishable by a fine dependent on the status of the victim. Ten-and-a-half cows per king murdered (no, really!), five for a major nobleman, and so on. At one particular event, we know the total fine levied, but not the number of victims. Chris wanted to find out the range of possibilities. This is an example of the change-making problem: given small change of various denominations, how can you give your customer 73p? I most often encounter this problem with rugby scores: if Wales beat England 67 - 3 (I can dream...), then at seven points for a converted try, five for an unconverted try and three for a penalty, what might have happened?

Chris had already written a Perl program that looked very much like this:
#!/usr/bin/perl                                                                 

use warnings;
use strict;

my $target = 23; # or whatever it really was
my $w_val = 10;
my $x_val = 5;
my $y_val = 2;
my $z_val = 1;

for my $w (0 .. $target / $w_val) {
    for my $x (0 .. $target/$x_val) {
        for my $y (0 .. $target / $y_val) {
            for my $z (0 .. $target / $z_val) {
                if ($w*$w_val + $x*$x_val + $y*$y_val + $z*$z_val == $target) {
                    print "$w, $x, $y, $z\n";
                }
            }
        }
    }
}

There are a few things to criticise in this program: it searches a lot of possibilities it can rule out (once $w and $x are too big, it doesn't matter what $y and $z are), it repeats a lot of calculations in the inner loops (more on that later), and the number of different social classes is baked into the structure of the program - if the early Irish lawmakers ever decide to institute a new class of noble, we'd have to modify our code rather than just changing the values of some parameters. But it's simple, easy to read, obviously correct and easy to write.

Chris wanted to re-write this in OCaml, which is a (mostly) functional language, and iteration is considered poor style by functional devotees. Could I help him re-write it to use recursion instead?

At this point I'm afraid I seized on the word "recursion", because, you see, there's an elegant recursive solution to the change-making problem and I thought he was asking about that. Somewhere between my wanting to entirely re-write the program using a new algorithm, it being trickier to write than I'd expected, and my not actually knowing OCaml very well, I managed to completely confuse him. Sorry, Chris. Anyway, we'll come back to the first program we wrote in a bit; here's what I think he was actually asking for:
let rec range i j = if i > j then [] else i :: (range (i+1) j)

let poss_range target value = range 0 (target/value)

let concat_map f xs = List.concat (List.map f xs)

let make_quad w x y z = (w,x,y,z)

let is_good target quad =
  match quad with
  (w,x,y,z) -> (w * 10 + x * 5 + y * 2 + z * 1) = target

let ways target = List.filter (is_good target)
  (concat_map
    (fun w -> concat_map
      (fun x -> concat_map
        (fun y -> List.map
          (make_quad w x y)
          (poss_range target 1))
        (poss_range target 2))
      (poss_range target 5))
    (poss_range target 10))

I think that's pretty clear: the only slightly tricksy bit is the use of the partially-applied function make_quad w x y. It's a couple of lines shorter than the original Perl version, but suffers from all of its flaws. It would be fairly easy to reduce the size of the search space - change the first argument passed to poss_range - but it's more interesting to fix the restriction to a fixed set of denominations. Here's an attempt to do that:
(* General utilities *)

let rec range i j = if i > j then [] else i :: (range (i+1) j)

let cons x xs = x::xs

let singleton x = [x]

let sum = List.fold_left (fun x y -> x + y) 0


(* dot product - x1*y1 + x2*y2 + x3*y3 + ... *)
let dot xs ys = sum (List.map2 (fun x y -> x * y) xs ys)

(* finite choice functions - ways of choosing one elt from each list *)
let rec choices xss =
  match xss with
    [] -> []
  | [xs] -> List.map singleton xs
  | xs::rest -> List.concat (List.map (fun x -> List.map (cons x) (choices rest)) xs)

(* functions specific to the problem at hand *)

let ranges target values = List.map (fun x -> range 0 (target / x)) values

let candidates target values = choices (ranges target values)

let isgood target values xs = (dot values xs) = target

let ways target values = List.filter (isgood target values) (candidates target values)

It doesn't look much like the previous version, but it works in exactly the same way: we iterate over all possible combinations of victim numbers, and check each one (using List.filter) to see if it gives the desired result. The meat of the program is the function choices, which generates finite choice functions. Given a list of lists [xs1;xs2;xs3;xs4;...] it will produce all lists whose first element is chosen from xs1, whose second element is chosen from xs2, etc. I believe it's considered poor style to use naked recursion like that (as opposed to using a higher-order function like List.map or List.fold_left), but frankly that function was the hardest thing I had to write in this whole exercise, and I was happy just to get something that worked.

Now, let's look at a more efficient approach.

The key to solving this problem efficiently is to realise that choosing a "first move" leaves us with a smaller example of the same type of problem. Let's suppose, as above, that Wales scored 67 points, and let's suppose that we've heard that they scored four converted tries. To work out the full space of possibilities, we just need to calculate all the ways of making 67 - 4*7 = 39 points with unconverted tries and penalties. But that's also a change-making problem! This means we can solve the problem recursively, and that's what I thought Chris was asking about. If we only have one denomination left to play with, then there's at most one way to solve the problem: return the list containing only target/value if value divides target evenly or the empty list (to indicate "no solutions") if it doesn't. If we have more than one denomination, we make all possible choices of "first move", and then recursively solve all the new, smaller problems that each choice presents us with.

Here's the code that we wrote to do this:
(* things that are probably in a standard library, if only I knew where to look
for them *)
let rec range i j = if i > j then [] else i :: (range (i+1) j)

let cons x xs = x::xs

let fuzz = 0.00001

(* from http://www.podval.org/~sds/ocaml-sucks.html *)
let round x = int_of_float (floor (x +. 0.5))

let float_eq x y = abs_float (x -. y) <= fuzz

let divides n d = let q = n /. d in float_eq q (float_of_int (truncate q))

(* the actual meat of the code *)
let rec ways ps tgt =
  match ps with
    [] -> []
  | [p] -> if (divides tgt p && tgt >= 0.) then [[round (tgt /. p)]] else []
  | p::rest -> List.concat (List.map (helper p rest tgt)
                  (range 0 (round (tgt /.p))))
and
  helper p ps tgt n =
    let ws = ways ps (tgt -. (float_of_int n) *. p) in
      List.map (cons n) ws

(* test code *)
let values = [10.;5.;2.;1.]

let target = 23.

let answers = ways values target

let add x y = x +. y

let dot xs ys = List.fold_left add 0.
  (List.map2 (fun x y -> x *. (float_of_int y)) xs ys)

let sums = List.map (dot values) answers

let errors = List.length (List.find_all (fun x -> not (float_eq x target)) sums)

let _ = print_endline ((string_of_int errors) ^ " errors.")

[One thing I didn't mention - though I can't remember the exact values in Chris' original problem, the death-price for a king really was 10.5 cows. Hence our program had to handle fractional denominations, which makes it a lot uglier than it would otherwise need to be - OCaml's lack of polymorphism means that it needs separate operators (called +., *. etc.) for floating-point numbers, and we need to strew conversion functions about the place.]

A refinement of this technique is to memoize the recursive function - this typically helps a lot, because you need to solve many of the same inner subproblems again and again. If you use the memoized version, the technique's called dynamic programming (a name chosen purely because it sounded impressive), which is a useful technique in all kinds of search applications. Some quick experiments indicate that it's not a win in this case, though - the memory requirements grow too rapidly.

Could we have written the recursive version in Perl? Sure:
#!/usr/bin/perl

use strict;
use warnings;

sub ways {
    my $target = shift;
    my $first = shift;
    if (scalar(@_) == 0 && $target >= 0 && ($target % $first == 0)) {
        return ([$target/$first]);
    }
    my @ways;
    for my $move (0 .. ($target/$first)) {
        push @ways, [$move, @$_] foreach ways($target - $move * $first, @_);
    }
    return @ways;
}
foreach (ways(23, 10,5,2,1)) { print join ",", @$_)."\n" };

[This version, by the way, was easily the most fun to write.]

I think it's interesting how different the equivalent Perl and OCaml versions look in this case: they're about the same length, but the Perl version is one subroutine, and the Caml version's divided into several smaller routines. This is because I find long subroutines in functional languages to be devilishly hard to write correctly - splitting expressions out into their own top-level functions makes it easier to ask the compiler what types it's inferring for them.

All suggestions for improvement, as ever, gratefully accepted :-)
Thursday, August 26th, 2010 03:41 pm (UTC)
If we only have one denomination left to play with, then there's at most one way to solve the problem

This is an optimization: for clarity, you might omit it, and let "no denominations left" be your own only base case.

I don't have OCaml lying around, but here's a first pass at the problem in Lisp, which should be similar, if slightly more verbose. Note also how I've avoided the need to do a "map" by building up a partial solution in the arguments to the recursive calls.
(defun ways (amount coins &optional prepend (count 0))
  "Return all ways of making change for AMOUNT using COINS.
AMOUNT is an integer and COINS a list of integers.
The list PREPEND is prepended to all solutions.
COUNT is added to the number of times the first coin is used."
  (if (endp coins)
      (if (= amount 0) (list prepend) nil)
    (let* ((coin (first coins))
           (unused (ways amount (rest coins) (append prepend (list count)))))
      (if (< amount coin)
          unused
        (append unused (ways (- amount coin) coins prepend (1+ count)))))))

(ways 23 '(10 5 2 1))
  ==> ((0 0 0 23) (0 0 1 21) (0 0 2 19) (0 0 3 17) (0 0 4 15) (0 0 5 13) (0 0 6 11) ...)
To an experienced Lisper, the presence of append in there is a bit of a red flag that there may be unnecessary allocation and list scanning going on. It's straightforward to eliminate one of the calls to append by putting even more of the solution into the arguments:
(defun ways (amount coins &optional prepend (count 0) more)
  "Return all ways of making change for AMOUNT using COINS.
AMOUNT is an integer and COINS a list of integers.
The list PREPEND is prepended to solutions.
COUNT is added to the number of times the first coin is used.
MORE is appended to the list of solutions."
  (if (endp coins)
      (if (= amount 0) (cons prepend more) more)
    (let* ((coin (first coins)))
      (ways amount (rest coins) (append prepend (list count)) 0
            (if (< amount coin)
                more
              (ways (- amount coin) coins prepend (1+ count) more))))))
This also has the benefit of putting one of the recursive calls into tail position, which will help the optimizer. Your mileage may differ as to whether the increase in speed is worth the decrease in readability.

We could eliminate the other call to append if we were willing to get reversed solutions (i.e. with the counts in opposite order to the coins list). I guess this depends on what the output is going to be used for.
Thursday, August 26th, 2010 04:22 pm (UTC)
Adding auxiliary arguments is a pretty standard technique in functional programming. The classic example is the factorial function (here, in Standard ML):
fun factorial 0 = 1
  | factorial n = n * factorial (n - 1)
which can be made tail-recursive by transforming it to
fun factorial n = 
    let fun factorial* 0 x = x
          | factorial* n x = factorial* (n-1) (n * x)
    in factorial* n 1
    end
Incidentally, when you have a fixed set of coins, and you're only interested in the number of ways to make change, there's a mathematical closed form for the answer. See, for example, Concrete Mathematics pages 344–346 ("the number of ways to change $1,000,000 is ... 66666793333412666685000001").
Thursday, August 26th, 2010 04:45 pm (UTC)
I'd encountered the general technique, but your use of it with count is IMHO rather elegant.

Thanks for the reference - I'll try to dig up a copy.
Thursday, August 26th, 2010 06:48 pm (UTC)
Ooh, thank you for this post, [livejournal.com profile] gareth_rees! I'm trying to learn to Lisp at the moment, and your code's really helpful!