In the first part of the series, I examined Newton-Lagrange interpolation, using functional language idioms. So far only F# have been put to the pit, but now Ruby and C# have dropped in. Geek alert: this post might become text-intensive and boring, however looking at the code itself should be enough to get the groove.
Ruby
This is everyone's favorite toy. And it has a good reason to be - it has most of the idioms functional languages use (and you'll see that by comparing the F# code), it has a lot of syntactic sugar, advanced constructs like meta-programming and closures, and all of that with minimal requirements from the programmer, in short - it stands out of your way very effectively. It flows. The first time I've tried Ruby was with the first Ruby on Rails rush (several years ago).
So, lets jump to the code and right after it I'll briefly explain things.
include Math
def newton_lagrange(xdata, ydata, px)
divdif = proc do | s, e |
if s == e
ydata[s]
else
(divdif.call(s+1, e) - divdif.call(s, e-1)) /
(xdata[e]-xdata[s])
end
end
divdifs = (0 ... xdata.length).to_a.map{ |x| divdif.call(0,x)}
xdata.reverse.zip(divdifs.reverse).inject(0){|res, (x, div)| res*(px-x)+div}
end
xdata = []; (-5).upto(5) { |x| xdata << x };
ydata = xdata.map { |x| 1.0 / (cos(2.0*x) + x**2.0)}
p newton_lagrange(xdata, ydata, 3.16)
So there are several constructs here, that compensate and fill up for F#'s capabilities with functional programming. Firstly and most importantly are Ruby Blocks. In ruby (and Scala, as we will see later) you can attach closures to code sections.
This means, that the code, internally, can access its attached codeblock and execute it at will. A codeblock can also define an anonymous function/lambda/closure for you, which leads us to the first section of the code.
divdif = proc do | <params> |
....<codeblock>
end
divdif is our familiar divided differences procedure from earlier on. Over here we're defining it as a closure, and we're explicitly calling it inside, since this is how we would imitate a letrec.
The second portion of the code is nifty and easy.
We activate our closure by mapping 0..K where K is some integer of our range, then we'll take each value and execute divdiv(0, val) just like we did in F# by currying.
Next comes the part thats supposed to compensate for F#'s advanced folding abilities.
We essentially take two lists, zip ((1 2), (3 4) => ((1 3) (2 4)) them together, and do a regular one-list fold (called inject in Ruby).
Last but not least is Ruby's range-generation syntax: num.upto(othernum) which is natural and cool. One more note, in the code we can see another syntactical sugar, {} attached to the same line, is the same codeblock we talked about just presented differently.
C# (3.0)
I've put 3.0 in parens, because even if I did or did not use its features, I belive its still a considerate step towards functionalism that it deserves its own respect.
So, lets take a look at the code, this is more verbose, but we'll discuss it after the break. Another note, beware of usage of Extension methods which may not seem obvious at first.
1 class Program
2 {
3 delegate double NumericDelegate(int s, int e);
4 static void Main(string[] args)
5 {
6 List<double> Y = new List<double>();
7 List<double> X;
8
9 //see also Sequence.Range.
10 X = (-5.0).UpTo(5.0);
11 X.ForEach(x => Y.Add(1.0 / (Math.Cos(2.0 * x) + Math.Pow(x, 2.0))));
12
13
14 Console.WriteLine(NewtonLagrange(X, Y, 3.16));
15 }
16
17 private static double NewtonLagrange(List<double> X, List<double> Y, double px)
18 {
19 List<double> divdifs = Enumerable.Range(0, X.Count).ToList().Select(x => divdif(0, x, X, Y)).ToList();
20
21 divdifs.Reverse();
22 X.Reverse();
23 IEnumerable<double[]> zipped = X.Zip(divdifs);
24 return zipped.Aggregate(0.0, (res, tup) => res * (px - tup[0]) + tup[1]);
25 }
26
27 private static double divdif(int s, int e, List<double> X, List<double> Y)
28 {
29 if (s == e)
30 return Y[s];
31 else
32 return (divdif(s + 1, e, X, Y) - divdif(s, e - 1, X, Y)) / (X[e] - X[s]);
33 }
34
35 }
36
37
38 public static class Exts
39 {
40 public static List<TSource []> Zip<TSource>(this IEnumerable<TSource> first, IEnumerable<TSource> second)
41 {
42 List<TSource[]> result = new List<TSource[]>();
43 IEnumerator<TSource> secondenum = second.GetEnumerator();
44 secondenum.MoveNext();
45 foreach (TSource item in first)
46 {
47 result.Add(new TSource[2] { item, secondenum.Current });
48 secondenum.MoveNext();
49 }
50
51 return result;
52 }
53
54 public static List<int> UpTo(this int x, int y, int step)
55 {
56 List<int> result = new List<int>();
57 while (x <= y)
58 {
59 result.Add(x);
60 x += step;
61 }
62 return result;
63 }
64 public static List<double> UpTo(this double x, double y, double step)
65 {
66 List<double> result = new List<double>();
67 while (x <= y)
68 {
69 result.Add(x);
70 x += step;
71 }
72 return result;
73 }
74 public static List<double> UpTo(this double x, double y)
75 {
76 return x.UpTo(y, 1.0);
77 }
78 }
First a note about the code. Some of the Extension method are redundant, and I've done it just to show how we could simulate Ruby's syntactic sugar.
It might take a bit longer to digest this code, because we're essentially moving from more declarative to less declarative (more imperative) programming style. So take the time to review.
Since the code is a bit long and more expressive, I'll make the explanation shorter .
19 List<double> divdifs = Enumerable.Range(0, X.Count).ToList().Select(x => divdif(0, x, X, Y)).ToList();
Here we're using Enumerable.Range, which was in fact Sequence.Range with earlier versions of the language. Doing a 'Select' is like doing a map. We can also see the new lambda syntax, <args> => { block }, I'm certain everyone is already fed up with talking about :).
24 return zipped.Aggregate(0.0, (res, tup) => res * (px - tup[0]) + tup[1]);
Aggregate is a 'fold' or 'inject'. To the keen eyed, you guessed right and this code is in fact almost identical to the Ruby code in structure and idioms, with missing parts implemented as needed, with the power of Extension Methods. I'll break this part here and let you all breath :).
Check back later as I'll put both of the code snippets up for download.
Next up: python and Scala.
Welcome to the Pit.
In this series I'll try and explore idiomatic functional language usages with C# 3.0, F#, Ruby (IronRuby), Python (IronPython), and Scheme. I will provide a short description for each snippet just to get motivated, and will close the series with performance and numerical stability analysis.
Test Subject
Like any good test, there has to be a good test subject. Now, while I could do a dull iteration / manipulation of a list of ice-cream flavors, I picked something a bit more interesting and practical: Newton-Lagrange interpolation.
Interpolation, in short, is taking partial data of something, and 'filling in the blanks' based on known context rules.
In elementary geometry, we all know that it takes 2 points in space to define a single unique line. But it also takes 3 points to define a single unique parabola -- a 2nd degree polynomial, 4 points for a 3rd degree polynomial, and so on.
The Newton-Lagrange method of interpolation will interpolate all of the missing data based on those rules, with ,of course, a certain error. In practice, Cubic Splines and Bezier curves (yea, thats the Pen tool in Photoshop), use forms of interpolation similar to that.
Theory
No, I'm not going to bore you with the math. Basically what needs to be done, is (preferably) compute a table of divided differences - here done with recursion, and use that with the Horner multiplication scheme.
F#
1 #light
2 open System
3
4 let NewtonLagrange (X, Y, px) =
5 let rec divdif s e =
6 match s, e with
7 | n1, n2 when n1 = n2 -> List.nth Y n1
8 | _, _ -> (divdif (s+1) e - divdif s (e-1) ) / (List.nth X e - List.nth X s)
9
10 let divdifs = List.map (divdif 0) [0 .. X.Length-1]
11 List.fold_right2 (fun a b acc -> acc*(px - a) + b) X divdifs 0.0
12
13
14
15 let X = [-5.0 .. 1.0 .. 5.0]
16 let Y = List.map (fun x -> 1.0 / (cos(2.0*x) + x**2.0)) X
17 NewtonLagrange(X, Y, 3.16)
let X = [-5.0 .. 1.0 .. 5.0]
Range: let X = [<start> .. <step> .. <end>]. Note, that in order to help F# do a good deal of inference, I like to state double values as 0.0 explicitly.
(fun x -> 1.0 / (cos(2.0*x) + x**2.0))
Lambdas: (fun <params> -> <body>), returns a function object. remember, functions are first-class citizens in functional langauges.
List.map (divdif 0) [0 .. X.Length-1]
Curried functions: This is a hyped, not that straight forward concept to grasp by the imperative programmer. A curried function is a function object, where we allow parameters to be partially specified. A curried function does not carry or curry anything, its name is a tribute to Haskell Curry. So, how could this be possible?
Imaging a function Add. Lets drop to a psuedo-F# syntax for a moment. So Add x y -> x+y, currying on it, we say AddFive = Add 5. AddFive now represent an add function that looks like: Add 5 y -> 5 + y. And thats all there is to it.
In line 10, we say: divdif 0. This goes and returns to map a divdif function where the parameter s is always zero.
List.fold_right2 (fun a b acc -> acc*(px - a) + b) X divdifs 0.0
List comprehensions: folding "folding" is just like "accumulating" or "injecting" or what ever idiom a language might use, to describe applying an operation to an item, and accumulating it, for every item on a list. As an example consider accumulating X = [ 1 1 1 1], which is folded into "4" by the "sum" function. Here we use a fold_right2, which starts from the end of the list, using 2 lists, which suits the Horner multiplication scheme perfectly.
let divdifs = List.map (divdif 0) [0 .. X.Length-1]
List comprehensions: map a map takes a list, applies a function to every item, and returns a 'deformed' list, that represents the same list with the function applied to every item.
An example would be taking a list X = [ 1 2 3 ] and squaring every element into [ 1 4 9 ].
Let rec, line 5,
a "let rec" is a binding operation, that will let you use the bound entity, inside its own body. In other words, its the way we do recursion here, and you'll also recognize that from Scheme, and other functional languages.
match s, e with
| decomposition1 when guard -> (handle case)
| decomposition2 -> (handle case)
Pattern matching: this feature would seem just like a sophisticated switch-statement, but its not, its a totally, completly, different concept. To grasp it fully, you might want to look at Prolog, or more practically, Erlang. In essence, a match block in F# can help you pick apart a matchable argument and do case-handling effectively, most preferably expressing recursion.
Paying a price: using List.nth is O(n)!. We can always use an Array and say myarray.(0), but I wanted to go pure functional anyway.
Next up: C#3.0 and Ruby.