Day to day stuff: October 2010

package nl.grons.whyfp
import scala.math.abs
/**

A translation of Miranda examples in chapter 4 of the paper Why Functional Programming Matters by John Hughes, 1984.

This article is the continuation of part 1 - Chapter 3 examples. I will show you lazy evaluation and how this can be useful.

@author Erik van Oosten
*/

object Chapter4_LazyList {

Chapter 4’s examples make heavy use of Miranda’s lazy evaluation. In particular, it uses the list as defined in chapter 3 in a lazy way. This allows the creation of lists of unbounded length. Scala supports lazy lists as well. You can find an implementation in the Stream class. However, as I want to stay close to the original, I redefine the list from chapter 3 as follows:

/** The list from chapter 3, redefined to make it a lazy list. */
sealed abstract class List[+A] {
  def head: A
  def tail: List[A]
}

/** Lazyly evaluate the argument of the constructor AND the value of tail. */
final class Cons[A](hd: A, tl: => List[A]) extends List[A] {
  override val head: A = hd
  override lazy val tail: List[A] = tl
}

If you compare this to the List implementation in the previous article, you note that parameter ‘tl’ of ‘Cons’ is passed as an by-name parameter (because of the =>). This means that the given expression is only evaluated when the parameter is used. Too prevent the tail expression to be evaluated again and again, the ‘tail’ value is declared as a val. Now ‘tl’ is evaluated exactly once. It is also declared as lazy so that the evaluation of the expression is delayed until the first time it is needed instead of during construction time.

These 3 things together allow the list to behave lazily.

You might also notice that ‘Cons’ is no longer a case class (not possible with by-name parameters). Therefore I define the following:

object Cons {
  def apply[A](hd: A, tl: => List[A]) = new Cons(hd, tl)
  // Warning: using unapply leads to evaluation of the tail.
  def unapply[A](c: Cons[A]) = Some(c.head, c.tail)
}

‘Nil’ did not change:

/** The empty list. */
case object Nil extends List[Nothing] {
  override def head: Nothing = {
    throw new NoSuchElementException("head of empty list")
  }
  override def tail: List[Nothing] = {
    throw new NoSuchElementException("tail of empty list")
  }
}

Just for fun: here is an example of how the list can be used to construct an infinite list.

/** A list of all integer numbers starting at x and higher. */
// Miranda: nat x = cons x nat(x+1)
def nat(x: Int): List[Int] = Cons(x, nat(x + 1))

The Miranda version of ‘nat’ is limited by the available memory, the Scala version is limited to scala.Int.MaxValue.

Function ‘map’ from the previous chapter won’t work. It will try to completely iterate the list which will obviously fail for unbounded lists. Lets implement it again with some simple pattern matching:

/** Create a new lazy list where each element a is replaced by f(a). */
def map_wrong[A, B](f: A => B)(l: => List[A]): List[B] = l match {
  case Nil => Nil
  case Cons(head, tail) => Cons(f(head), map(f)(tail))
}

This implementation seemed to work for quite some time, even though everything was slow. If you look carefully, you will see that it evaluates ‘tail’ in ‘Cons.unapply’ which is invoked to match the pattern Cons(head, tail). The result is that exactly one more item in the list is evaluated then necessary. Though this is merely wasteful for most algorithms, it is catastrophic (as in stack overflow) when that evaluation contains another recursive call to ‘map’. This actually happens in the last example of this article.

Here is an implementation that does not have this problem:

/** Create a new lazy list where each element a is replaced by f(a). */
def map[A, B](f: A => B)(l: => List[A]): List[B] = l match {
  case Nil => Nil
  case c: Cons[A] => Cons(f(c.head), map(f)(c.tail))
}

Notice that the recursive expression, which includes c.tail, is now used in a place that is lazily evaluated.

On the console:

scala> import nl.grons.whyfp.Chapter4_LazyList._

scala> def printFirst[A](count: Int, l: List[A]) {
  if (count > 0 && l != Nil) { println(l.head); printFirst(count - 1, l.tail) } }
printFirst: [A](count: Int,l: nl.grons.whyfp.Chapter4_LazyList.List[A])Unit

scala> printFirst(3, map[Int,Int](_ * 2)(nat(3)))
6
8
10

*/
}

object Chapter4_1_NewtonRaphsonSquareRoots {

import Chapter4_LazyList._
/*

Here’s an imperative Scala implementation that does a numerical calculation of the square root of a number.

/** Approximate the square root of n to within eps
 * starting with approximation a0. */
def imperative_squareroot(a0: Double, eps: Double, n: Double): Double = {
  // Warning: imperative scala!
  var x = a0
  // The initial value of y does not matter so long as abs(x-y) > eps
  var y = a0 + 2 * eps
  while (abs(x - y) > eps) {
    y = x
    x = (x + n/x) / 2
  }
  x
}

Here are my translations of the rest:

/** Next approximation of the square root of n,
 * based on the previous approximation x. */
// Miranda: next N x = (x + N/x) / 2
def next(n: Double)(x: Double): Double = (x + n/x) / 2

/** Produce an infinite list that starts with a and where all other values
 * are the result of applying f to the previous value. */
// Miranda: repeat f a = cons a (repeat f (f a))
def repeat[A](f: A => A, a: A): List[A] = Cons(a, repeat(f, f(a)))

/** Give the first value of the list that is within eps of its preceding value. */
// Miranda:
// within eps (cons a (cons b rest)) =
//         = b,                        if abs(a-b) <=  eps
//         = within eps (cons b rest), otherwise
def within(eps: Double, l: List[Double]): Double = {
  val a = l.head
  if (a.isNaN) throw new RuntimeException("nan")
  val rest = l.tail
  val b = rest.head
  if(abs(a - b) <= eps) b else within(eps, rest)
}

When a calculation overflow or underflow occurs (not unlikely when working with approximations), ‘within’ goes amok and will loop forever. This is prevented by validating the head of the list against NaN, Not A Number.

/** Approximate the square root of n to within eps
 * starting with approximation a0. */
// Miranda: sqrt a0 eps N = within eps (repeat (next N) a0)
def sqrt(a0: Double, eps: Double, n: Double) = within(eps, repeat(next(n), a0))

/** Gives the first value of the list that of which the ratio of change
 * with its preceding value is lower then eps. */
// Miranda:
// relative eps (cons a (cons b rest)) =
//         = b,                          if abs(a-b) <= eps*abs b
//         = relative eps (cons b rest), otherwise
def relative(eps: Double, l: List[Double]): Double = {
  val a = l.head
  if (a.isNaN) throw new RuntimeException("nan")
  val rest = l.tail
  val b = rest.head
  if(abs(a - b) <= eps * abs(b)) b else relative(eps, rest)
}

/** Approximate the square root of n to within ratio eps starting with
 * approximation a0. */
// Miranda: relativesqrt a0 eps N = relative eps (repeat (next N) a0)
def relativesqrt(a0: Double, eps: Double, n: Double) =
  relative(eps, repeat(next(n), a0))

On the Console:

scala> import nl.grons.whyfp.Chapter4_LazyList._
scala> import nl.grons.whyfp.Chapter4_1_NewtonRaphsonSquareRoots._

scala> relativesqrt(0.1, 0.01, 4)
res0: Double = 2.000010925778043

*/
}

object Chapter4_2_NumericalDifferentiation {

import Chapter4_LazyList._
import Chapter4_1_NewtonRaphsonSquareRoots._
/*

This chapter is about numerical differentiation.

/** An approximation of the differentiation of f for x over h. */
// Miranda: easydiff f x h = (f(x+h)-f x) / h
def easydiff(f: Double => Double, x: Double)(h: Double): Double =
  (f(x + h) - f(x)) / h

// Miranda:
// differentiate h0 f x = map (easydiff f x)(repeat halve h0)
// halve x = x/2
def differentiate(h0: Double, f: Double => Double, x: Double) =
  map(easydiff(f, x))(repeat(halve, h0))
def halve(x: Double) = x / 2

On the console:

scala> import nl.grons.whyfp.Chapter4_LazyList._
scala> import nl.grons.whyfp.Chapter4_1_NewtonRaphsonSquareRoots._
scala> import nl.grons.whyfp.Chapter4_2_NumericalDifferentiation._

scala> val f: Double => Double = math.pow(_, 3.0)
f: (Double) => Double = <function1>

scala> within(0.001, differentiate(1, f, 2))
res0: Double = 12.000732436776161

Some more examples:

// Miranda:
// elimerror n (cons a (cons b rest)) =
//        = cons ((b*(2**n)-a)/(2**n-1)) (elimerror n (cons b rest))
def elimerror(n: Double, l: List[Double]): List[Double] = {
  val a = l.head
  val bandrest = l.tail
  val b = bandrest.head
  Cons(
    (b * math.pow(2.0, n) - a) / (math.pow(2.0, n) - 1.0),
    elimerror(n, bandrest))
}

/** Calculate the order of n to be used in elimerror. */
// Miranda:
// order (cons a (cons b (cons c rest))) =
//         = round(log2( (a-c)/(b-c) - 1 ))
// round x = x rounded to the nearest integer
// log2 x = the logarithm of x to the base 2
def order(l: List[Double]): Double = {
  val LOG2: Double = math.log(2.0)
  def log2(x: Double) = math.log(x) / LOG2

  val a = l.head
  val b = l.tail.head
  val c = l.tail.tail.head
  val o = math.round(log2((a-c)/(b-c) - 1))
  if (o.isNaN || o == 0) 1 else o
}

// Miranda: improve s = elimerror (order s) s
def improve(s: List[Double]): List[Double] = elimerror(order(s), s)

Note that my implementation of ‘order’ is not only real (no vague descriptions of ‘round’ and ‘log2’), but it is more robust then the original. First, it protects against NaN for the case that b-c reaches 0. Secondly, it prevents a result of 0 as that leads to a division by 0 in ‘elimerror’.

On the console:

scala> within(0.001, improve(differentiate(1, f, 2)))
res0: Double = 11.9998779296875

scala> // Improve can be improved recursively

scala> within(0.001, improve(improve(improve(differentiate(1, f, 2)))))
res1: Double = 12.0

// Miranda:
// super s = map second (repeat improve s)
// second (cons a (cons b rest)) = b
def superimprove(s: List[Double]): List[Double] = map(second)(repeat(improve, s))
def second(l: List[Double]): Double = l.tail.head

def superdifferentiate(f: Double => Double, x: Double): Double = {
  val eps: Double = 0.00000001
  within(eps, superimprove(differentiate(x + 2 * eps, f, x)))
}

}

object Chapter4_3_NumericalIntegration {

import Chapter4_LazyList._
/*

This chapter is about numerical integration.

// Miranda: easyintegrate f a b = (f a + f b)*(b-a)/2
def easyintegrate(f: Double => Double, a: Double, b: Double) =
  (f(a) + f(b)) * (b-a) / 2

// Miranda: integrate f a b = cons (easyintegrate f a b)
//                                 (map addpair (zip (integrate f a mid)
//                                                   (integrate f mid b)))
//                  where mid = (a+b)/2
def integrate_simple(f: Double => Double, a: Double, b: Double): List[Double] = {
  def addpair(pair: (Double, Double)) = pair._1 + pair._2
  val mid = (a + b) / 2
  Cons(
    easyintegrate(f, a, b),
    map(addpair)(zip(integrate(f, a, mid), integrate(f, mid, b)))
  )
}

/** Convert two lists to a list of pairs. */
// Miranda: zip (cons a s) (cons b t) = cons (pair a b) (zip s t)
def zip[A,B](as: => List[A], bs: => List[B]): List[(A,B)] = as match {
  case Nil => Nil
  case a: Cons[A] => bs match {
    case Nil => Nil
    case b: Cons[B] => Cons((a.head, b.head), zip(a.tail, b.tail))
  }
}

Again note that ‘zip’ was defined such that a.tail and b.tail are defined as a lazily evaluated argument.

// Miranda:
// integrate f a b = integ f a b (f a) (f b)
// integ f a b fa fb = cons ((fa+fb)*(b-a)/2)
//                          (map addpair (zip (integ f a m fa fm)
//                                            (integ f m b fm fb)))
//         where m = (a+b)/2
//               fm = f m
def integrate(f: Double => Double, a: Double, b: Double): List[Double] = {
  def integ(a: Double, b: Double, fa: Double, fb: Double): List[Double] = {
    def addpair(pair: (Double, Double)) = pair._1 + pair._2

    val mid = (a + b) / 2
    val fm = f(mid)
    Cons(
      (fa + fb) * (b-a) / 2,
      map(addpair)(zip(integ(a, mid, fa, fm), integ(mid, b, fm, fb)))
    )
  }

  integ(a, b, f(a), f(b))
}

And of course, seeing is believing:

scala> import nl.grons.whyfp.Chapter4_LazyList._
scala> import nl.grons.whyfp.Chapter4_1_NewtonRaphsonSquareRoots._
scala> import nl.grons.whyfp.Chapter4_2_NumericalDifferentiation._
scala> import nl.grons.whyfp.Chapter4_3_NumericalIntegration._

scala> def f(x: Double): Double = 1/(1+x*x)
f: (x: Double)Double

scala> within(0.001, integrate(f, 1, 2))
res0: Double = 0.3218612363325556

scala> relative(0.001, integrate(f, 1, 2))
res1: Double = 0.3217782239721789

scala> relative(0.001, superimprove(integrate(f, 1, 2)))
res2: Double = 0.3217525935931843

*/
}
/*

Conclusion

Scala is an excellent language for functional style programs and I really like that static typing prevents tons of errors. My only problem with Scala — in terms of the ‘Why Functional Programming matters’ paper — is that lazy evaluation is not the natural way of working. The way expressions are evaluation is more closely to the Java world. To make use of lazy evaluation you have to be very much aware of what your doing and carefully make use of the lazy keyword and by-name parameters. With Miranda this comes naturally though probably at the cost of some performance for the case laziness is not needed. Perhaps that the standard collections library (scala.collection.immutable.Stream) hides this sufficiently and makes life easy.

A more fundamental way to fix this problem would be to express the laziness nature of an expression in its type. I am out of the academic world for a looong time, so I have no idea if this is being researched at all. It will take some hard thoughts to make this practical and simple to use though.

The next challenge would be to repeat this experiment with the standard collections. But I’ll leave that to someone with more time...

Day to day stuff

Tuesday, October 5, 2010

Why Functional Programming Matters - An exploration of functional Scala - Part 2

object Chapter4_LazyList {

object Chapter4_1_NewtonRaphsonSquareRoots {

object Chapter4_2_NumericalDifferentiation {

object Chapter4_3_NumericalIntegration {

Conclusion