Functional-programming.rmd

---
title: Functional programming
layout: default
---

# Functional programming

At its heart, R is a functional programming (FP) language; it focusses on the creation and manipulation of functions. R has what's known as first class functions, functions that can be:

* created without a name,
* assigned to variables and stored in lists,
* returned from functions, and
* passed as arguments to other functions.

This means that you can do anything with functions that you can do with vectors: you can create them inside other functions, pass them as arguments to functions, return them as results from functions and store multiple functions in a list. This chapter will explore the consequences of R's functional nature and introduce a new set of techniques for removing redundancy and duplication in your code. We'll start with a motivating example, showing how you can use functional programming techniques to reduce duplication in some typical code for cleaning data and summarising data. This example will introduce some of the key building blocks of functional programming, which we will then dive into in more detail:

* __Anonymous functions__, functions that don't have a name

* __Closures__, functions written by other functions

* __Lists of functions__, storing functions in a list

The chapter concludes with a case study exploring __numerical integration__ showing how we can build a family of composite integration tools starting from very simple primitives.  This will be a recurring theme: if we start with small building blocks that we can easily understand, when we combine them into more complex structures, we can still feel confident that they are correct.

The exposition of functional programming continues in the following two chapters: [[functionals]], which explore functions that take functions as arguments and give vectors as output, and [[function operators]], functions that both input and output functions.

### Other languages

FP programming techniques are the core technique in FP languages, like Haskell, OCaml and F#. They are also well supported in multi-paradigm systems like Lisp, Scheme, Clojure and Scala. You can use FP techniques in modern scripting languages, like python, ruby and javascript, but they tend not to be the dominant technique employed by most programmers. Java and C# provide few functional tools, and while it's possible to do FP in those languages, it tends to be a somewhat awkward fit. Similarly, for functional programming in C. Googling for "functional programming in X" will find you a tutorial in any language, but it may be syntactically awkward or used so rarely that other programmers will not understand your code.

Recently FP has experienced a surge in interest because it provides a complementary set of techniques to object oriented programming, which has been the dominant style for the last several decades. Since FP functions tend to not modify their inputs, it makes programs that are easier to reason about using only local information, and are often easier to parallelise. The traditional weaknesses of FP languages, poorer performance and sometimes unpredictable memory usage, have been largely eliminated in recent years.

## Motivation

Imagine you've loaded a data file that uses -99 to represent missing values, like the following sample dataset.

```R
# Generate a sample dataset
set.seed(1014)
df <- data.frame(replicate(6, sample(c(1:10, -99), 10, rep = T)))
names(df) <- letters[1:6]
head(df)
```

When you first start writing R code, you might write code like the following, dealing with the duplicated processing of each column with copy-and-paste:

```R
df$a[df$a == -99] <- NA
df$b[df$b == -99] <- NA
df$c[df$c == -98] <- NA
df$d[df$d == -99] <- NA
df$e[df$e == -99] <- NA
df$f[df$g == -99] <- NA
```

One problem with copy-and-paste is that it's easy to make mistakes: can you spot two in the block above? The problem is that one idea, that missing values are represent as -99, is repeated many times. Repitition is bad because it allows for inconsistencies (aka bugs), and it makes the code harder to change. For example, if the representation of missing values changes from -99 to 9999, then we need to make the change in many places, not just one.

The "do not repeat yourself", or DRY, principle, was popularised by the [pragmatic programmers](http://pragprog.com/about), Dave Thomas and Andy Hunt. This principle states that "every piece of knowledge must have a single, unambiguous, authoritative representation within a system". Adhering to this principle prevents bugs caused by inconsistencies, and makes software that is easier to adapt to changing requirements. The ideas of FP are important because they give us new tools to reduce duplication.

We can start applying some of the ideas of FP to our example by writing a function that fixes the missing values in a single vector:

```R
fix_missing <- function(x) {
  x[x == -99] <- NA
  x
}
df$a <- fix_missing(df$a)
df$b <- fix_missing(df$b)
df$c <- fix_missing(df$c)
df$d <- fix_missing(df$d)
df$e <- fix_missing(df$e)
df$f <- fix_missing(df$e)
```

This reduces the scope for errors, but doesn't eliminate them.  We've still made an error, because we've repeatedly applied our function to each column. To prevent that error from occuring we need to remove the copy-and-paste application of our function to each column. To do this, we need to combine, or compose, our function for correcting missing values with a function that does something to each column in a data frame, like `lapply()`.  

`lapply()` takes three inputs: `x`, a list; `f`, a function; and ``...`, other arguments to pass to `f`. It applies the function to each element of the list and returns a new list. Since data frames are also lists, `lapply()` also works on data frames. `lapply(x, f, ...)` is equivalent to the following for loop:

```R
out <- vector("list", length(x))
for (i in seq_along(x)) {
  out[[i]] <- f(x[[i]], ...)
}
```

The real `lapply()` is rather more complicated since it's implemented in C for efficiency, but the essence of the algorithm is the same. `lapply()` is called a __functional__, because it takes a function as an arguments. Functionals are an important part of functional programming and we'll learn more about them in the [[functionals]] chapter.

We can use `lapply()` with one small trick: rather than simply assigning the results to `df` we assign them to `df[]`, so R's usual subsetting rules take over and we get a data frame instead of a list. (If this comes as a surprise, you might want to read over the [[subsetting]] appendix)

```R
fix_missing <- function(x) {
  x[x == -99] <- NA
  x
}
df[] <- lapply(df, fix_missing)
```

As well as being more compact, there are four main advantages of this code over our previous code:

* If the representation of missing values changes, we only need to change it in one place.

* There is no way for some columns to be treated differently than others.

* Our code works regardless of the number of columns in the data frame, and there is no way to miss a column because of a copy and paste error.

* It is easy to generalise this technique to a subset of our columns:

    ```R
    df[1:5] <- lapply(df[1:5], fix_missing)
    ```

The key idea here is composition. We take two simple functions, one which does something to each column, and one which fixes missing values, and combine them together to fix missing values in every column. Writing simple functions than can be understood in isolation and then composed together to solve complex problems is an important technique for effective FP.

What if different columns used different indicators for missing values? You again might be tempted to copy-and-paste:

```R
fix_missing_99 <- function(x) {
  x[x == -99] <- NA
  x
}
fix_missing_999 <- function(x) {
  x[x == -999] <- NA
  x
}
fix_missing_9999 <- function(x) {
  x[x == -999] <- NA
  x
}
```

But as before, it's easy to create bugs. The next functional programming tool we'll talk about helps deal with this sort of duplication: when we have multiple functions that all follow the same basic template. Closures, functions that return functions, allow us to make many functions from a template:

```R
missing_fixer <- function(na_value) {
  function(x) {
    x[x == na_value] <- NA
    x
  }
}
fix_missing_99 <- missing_fixer(-99)
fix_missing_999 <- missing_fixer(-999)
fix_missing_9999 <- missing_fixer(-9999)
```

(In this case, you could argue that we should just add another argument:

```R
fix_missing <- function(x, na.value) {
  x[x == na.value] <- NA
  x
} 
```

That's a reasonable solution here, but it doesn't work so well in every situation. We'll see more compelling uses for closures later in the chapter.)

Let's now consider a new problem: once we've cleaned up our data, we might want to compute the same set of numerical summaries for each variable. We could write code like this:

```R
mean(df$a)
median(df$a)
sd(df$a)
mad(df$a)
IQR(df$a)

mean(df$b)
median(df$b)
sd(df$b)
mad(df$b)
IQR(df$b)
```

But we'd be better off identifying the sources of duplication and then removing them. Take a minute or two to think about how you might tackle this problem before reading on.

One approach would be to write a summary function and then apply it to each column:

```R
summary <- function(x) { 
  c(mean(x), median(x), sd(x), mad(x), IQR(x))
}
lapply(df, summary)
```

But there's still some duplication here. If we make the summary function slightly more realistic, it's easier to see the duplication:

```R
summary <- function(x) { 
 c(mean(x, na.rm = TRUE), 
   median(x, na.rm = TRUE), 
   sd(x, na.rm = TRUE), 
   mad(x, na.rm = TRUE), 
   IQR(x, na.rm = TRUE))
}
```

All five functions are called with the same arguments (`x` and `na.rm`) which we had to repeat five times. As before, this duplication makes our code fragile: it makes it easier to introduce bugs and harder to adapt to changing requirements. 

We can take advantage of another functional programming technique, storing functions in lists, to remove this duplication:

```R
summary <- function(x) {
  funs <- c(mean, median, sd, mad, IQR)
  lapply(funs, function(f) f(x, na.rm = TRUE))
}
```

The remainder of this chapter will discuss these techniques in more detail. But before we can start on those more complicated techniques, we need to start by revising a simple functional programming tool, anonymous functions.

## Anonymous functions

In R, functions are objects in their own right. They aren't automatically bound to a name and R doesn't have a special syntax for creating named functions, unlike C, C++, python or ruby. You might have noticed this already, because when you create a function, you use the usual assignment operator to give it a name. 

Given the name of a function, like `"mean"`, it's possible to find the function using `match.fun()`. You can't do the opposite: given the object `f <- mean`, there's no way to find its name. Not all functions have a name, and some functions have more than one name. Functions that don't have a name are called __anonymous functions__.

We use anonymous functions when it's not worth the effort of creating a named function:

```R
lapply(mtcars, function(x) length(unique(x)))
Filter(function(x) !is.numeric(x), mtcars)
integrate(function(x) sin(x) ^ 2, 0, pi)
```

Unfortunately the default R syntax for anonymous functions is quite verbose. To make things a little more concise, the `pryr` packages provides `f()`:

```R
lapply(mtcars, f(length(unique(x))))
Filter(f(!is.numeric(x)), mtcars)
integrate(f(sin(x) ^ 2), 0, pi)
```

I'm not still sure whether I like this style or not, but it sure is compact!  Other similar ideas are implemented in `gsubfn::fn()` and `ptools::fun()`.

Like all functions in R, anoynmous functions have `formals()`, a `body()`, and a parent `environment()`:
  
```R
formals(function(x = 4) g(x) + h(x))
body(function(x = 4) g(x) + h(x))
environment(function(x = 4) g(x) + h(x))
```

You can call anonymous functions directly, but the code is a little tricky to read because you must use parentheses in two different ways: to call a function, and to make it clear that we want to call the anonymous function `function(x) 3`, not inside our anonymous function call a function called `3` (which isn't a valid function name!):

```R
(function(x) x + 3)(10)

# Exactly the same as
f <- function(x) x + 3
f(10)

# Doesn't do what you expect
function(x) 3()
```

You can supply arguments to anonymous functions in all the usual ways (by position, exact name and partial name) but if you find yourself doing this, it's probably a sign that your function needs a name.

One of the most common uses for anonymous functions is to create closures, functions made by other functions. Closures are described in the next section.

### Exercises

* Use `lapply()` and an anonymous function to find the coefficient of variation (the standard deviation divided by the mean) for all columns in the `mtcars` dataset

* Use `integrate()` and an anonymous function to find the area under the curve of:

  * `y = x ^ 2 - x`, x in [0, 10]
  * `y = sin(x) + cos(x)`, x in [-pi, pi]
  * `y = exp(x) / x`, x in [10, 20]
  
  Use [wolframalpha](http://www.wolframalpha.com/) to check your answers.

* A good rule of thumb is that an anonymous function should fit on one line and shouldn't need to use `{}`.  Review your code: where could you have used an anonymous function instead of a named function? Where should you have used a named function instead of an anonymous function?

## Introduction to closures

"An object is data with functions. A closure is a function with data." 
--- [John D Cook](http://twitter.com/JohnDCook/status/29670670701)

One use of anonymous functions is to create small functions that it's not worth naming; the other main use of anonymous functions is to create closures, functions written by functions. Closures are so called because they __enclose__ the environment of the parent function, and can access all variables in the parent. This is useful because it allows us to have two levels of parameters. One level of parameters (the parent) controls how the function works; the other level (the child) does the work. The following example shows how we can use this idea to generate a family of power functions. The parent function (`power()`) creates child functions (`square()` and `cube()`) that do the work.

```R
power <- function(exponent) {
  function(x) x ^ exponent
}

square <- power(2)
square(2)
square(4)

cube <- power(3)
cube(2)
cube(4)
```

In R, almost every function is a closure, because all functions remember the environment in which they are created, typically either the global environment, if it's a function that you've written, or a package environment, if it's a function that someone else has written. The only exception are primitive functions, which call to C directly.

When you print a closure, you don't see anything terribly useful:

```R
square
cube
```

That's because the function itself doesn't change; it's the enclosing environment, `environment(square)`, that's different. One way to see the contents of the environment is to convert it to a list:

```R
as.list(environment(square))
as.list(environment(cube))
```

Another way to see what's going on is to use `pryr::unenclose()`, which substitutes the variables defined in the enclosing environment into the original functon:

```R
library(pryr)
unenclose(square)
unenclose(cube)
```

Note that the parent environment of the closure is the environment created when the parent function is called:

```R
power <- function(exponent) {
  print(environment())
  function(x) x ^ exponent
}
zero <- power(0)
environment(zero)
```

This environment normally disappears once the function finishes executing, but because we return a function, the environment is captured and attached to the new function. Each time we re-run `power()` a new environment is created, so each function produced by power is independent.

Closures are useful for making function factories, and are one way to manage mutable state in R. 

### Function factories

We've already seen two example of function factories, `missing_fixer()` and `power()`. In both these cases using a function factory instead of a single function with multiple arguments has little, if any, benefit. Function factories are most useful when:

* the different levels are more complex, with multiple arguments and complicated bodies

* some work only needs to be done once, when the function is generated

INSERT USEFUL EXAMPLE HERE

We'll see another compelling use of function factories when we learn more about [[functionals]]; they are particularly well suited to maximum likelihood problems.

### Mutable state

Having variables at two levels makes it possible to maintain state across function invocations, because while the function environment is refreshed every time, its parent environment stays constant. The key to managing variables at different levels is the double arrow assignment operator (`<<-`). Unlike the usual single arrow assignment (`<-`) that always assigns in the current environment, the double arrow operator will keep looking up the chain of parent environments until it finds a matching name. ([[Environments]] has more details on how it works)

Together, a static parent environment and `<<-` make it possible to maintain state across function calls.  The following example shows a counter that records how many times a function has been called. Each time `new_counter` is run, it creates an environment, initialises the counter `i` in this environment, and then creates a new function.

```R
new_counter <- function() {
  i <- 0
  function() {
    i <<- i + 1
    i
  }
}
```

The new function is a closure, and its enclosing environment is the usually temporary environment created with `new_counter` is run.  When the closures `counter_one` and `counter_two` are run, each one modifies the counter in a different enclosing environment and so maintain different counts. 

```R
counter_one <- new_counter()
counter_two <- new_counter()

counter_one() # -> [1] 1
counter_one() # -> [1] 2
counter_two() # -> [1] 1
```

We can use our environment inspection tools to see what's going on here:

```R
as.list(environment(counter_one))
as.list(environment(counter_two))
```

The counters get around the "fresh start" limitation by not modifying variables in their local environment. Since the changes are made in the unchanging parent (or enclosing) environment, they are preserved across function calls.

What happens if we don't use a closure? What happens if we only use `<-` instead of `<<-`? Make predictions about what will happen if you replace `new_counter()` with each variant below, then run the code and check your predictions.

```R
i <- 0
new_counter2 <- function() {
  i <<- i + 1
  i
}
new_counter3 <- function() {
  i <- 0
  function() {
    i <- i + 1
    i
  }
}
```

Modifying values in a parent environment is an important technique because it is one way to generate "mutable state" in R. Mutable state is normally hard to achieve, because every time it looks like you're modifying an object, you're actually creating a copy and modifying that. That said, if you do need mutable objects, except in the simplest of cases, it's usually better to use the [[R5]] OO system. R5 objects are easier to document, and provide easier ways to inherit behaviour.

The power of closures is tightly coupled to [[functionals]] and [[function operators]], and you'll see many more examples of closures in those two chapters.  The following section disucsses the remaining important property of functions: the ability to store them in a list.

### Exercises

* What does the following statistical function do? What would be a better name for it? (The existing name is a bit of a hint)

    ```R
    bc <- function(lambda) {
      if (lambda == 0) {
        function(x) log(x)
      } else {
        function(x) (x ^ lambda - 1) / lambda
      }
    }
    ```

* Create a function that creates functions that compute the ith [central moment](http://en.wikipedia.org/wiki/Central_moment) of a numeric vector. You can test it by running the following code:

    ```R
    m1 <- moment(1)
    m2 <- moment(2)

    x <- runif(m1, 100)
    stopifnot(all.equal(m1(x), mean(x)))
    stopifnot(all.equal(m2(x), var(x) * 99 / 100))
    ```

* What does `approxfun()` do? What does it return? 

* What does the `ecdf()` function do? What does it return? 

* Create a function `pick()`, that takes an index, `i`, as an argument and returns a function an argument `x` that subsets `x` with `i`.
  
    ```R
    lapply(mtcars, pick(5))
    # should do the same this as
    lapply(mtcars, function(x) x[[5]])
    ```

## Lists of functions

In R, functions can be stored in lists. Instead of giving a set of functions related names, you can store them in a list.  This makes it easier to work with groups of related functions, in the same way a data frame makes it easier to work with groups of related vectors.

We'll start with a simple example: benchmarking, when you are comparing the performance of multiple approaches to the same problem. For example, if you wanted to compare a few approaches to computing the mean, you could store each approach (function) in a list:

```R
compute_mean <- list(
  base = function(x) mean(x),
  sum = function(x) sum(x) / length(x),
  manual = function(x) {
    total <- 0
    n <- length(x)
    for (i in seq_along(x)) {
      total <- total + x[i] / n
    }
    total
  }
)
```

Calling a function from a list is straightforward: just get it out of the list first:

```R
x <- runif(1e5)
system.time(compute_mean$base(x))
system.time(compute_mean[[2]](x))
system.time(compute_mean[["manual"]](x))
```

If we want to call each functions to check that we've implemented them correctly and they return the same answer, we can use `lapply()`, either with an anonymous function, or an equivalent named function.

```R
lapply(compute_mean, function(f, ...) f(...), x)

call_fun <- function(f, ...) f(...)
lapply(compute_mean, call_fun, x)
```

If we want to time how long each function takes, we can combine lapply with `system.time()`:

```R
lapply(compute_mean, function(f) system.time(f(x)))
```
    
Coming back to our original motivating example, another use case of lists of functions is summarising an object in multiple ways. We could store each summary function in a list, and then run them all with `lapply()`:

```R
funs <- list(
  sum = sum,
  mean = mean,
  median = median
)
lapply(funs, function(f) f(1:10))
```

What if we wanted our summary functions to automatically remove missing values?  One approach would be make a list of anonymous functions that call our summary functions with the appropriate arguments:

```R
funs2 <- list(
  sum = function(x, ...) sum(x, ..., na.rm = TRUE),
  mean = function(x, ...) mean(x, ..., na.rm = TRUE),
  median = function(x, ...) median(x, ..., na.rm = TRUE)
)
```

But this leads to a lot of duplication - each function is almost identical apart from a different function name. We could write a closure to abstract this away:

Instead, we could modify our original `lapply()` call:

```R
lapply(funs, function(f) f(x))
lapply(funs, function(f) f(x, na.rm = TRUE))

# Or use a named function instead of an anonymous function
remove_missings <- function(f) {
  function(...) f(..., na.rm = TRUE)
}
funs2 <- lapply(funs, remove_missings)
```

### Moving lists of functions to the global environment

From time to time you may want to create a list of functions that you want to be available to your users without having to use a special syntax. For a simple example, imagine you want to make it easy to create HTML with code, by mapping each html tag to an R function. The following simple example creates functions for `<p>` (paragraphics), `<b>` (bold), `<i>` (italics), and `<img>` (images).  Note the use of a closure function factory to produce the text for `<p>`, `<b>` and `<i>` tags.

```R
simple_tag <- function(tag) {
  function(...) paste0("<", tag, ">", paste0(...), "</", tag, ">")
}
html <- list(
  p = simple_tag("p"),
  b = simple_tag("b"),
  i = simple_tag("i"),
  img = function(path, width, height) {
    paste0("<img src='", path, "' width='", width, "' height = '", height, '" />')
  }
)
```

We store the functions in a list because we don't want them to be available all the time: the risk of a conflict between an existing R function and an HTML tag is high. However, keeping them in a list means that our code is more verbose than necessary:

```R
html$p("This is ", html$b("bold"), ", ", html$i("italic"), " and ",
   html$b(html$i("bold italic")), " text")
```

We have three options to eliminate the use of `html$`, depending on how long we want the effect to last:

* For a very temporary effect, we can use a `with()` block:

    ```R
    with(html, p("This is ", b("bold"), ", ", i("italic"), " and ",
      b(i("bold italic")), " text"))
    ```

* For a longer effect, we can use `attach()` to add the functions in `html` in to the search path.  It's possible to undo this action using `detach`:

    ```R
    attach(html)
    p("This is ", b("bold"), ", ", i("italic"), " and ",
      b(i("bold italic")), " text")
    detach(html)
    ```

* Finally, we could copy the functions into the global environment with `list2env()`.  We can undo this action by deleting the functions after we're done.
  
    ```R
    list2env(html, environment())
    p("This is ", b("bold"), ", ", i("italic"), " and ",
      b(i("bold italic")), " text")
    rm(list = names(html), envir = environment())
    ```

I recommend the first option because it makes it very clear what's going on, and when code is being executed in a special context.

### Exercises

* Implement a summary function that works like `base::summary()`, but takes a list of functions to use to compute the summary. Modify the function so it returns a closure, making it possible to use it as a function factory.

* Create a named list of all base functions.  Use `ls()`, `get()` and `is.function()`. Use that list of functions to answer the following questions:

    * Which base function has the most arguments? 
    * How many base functions have no arguments?

* Which of the following commands is `with(x, f(z))` equivalent to?

    (a) `x$f(x$z)`
    (b) `f(x$z)`
    (c) `x$f(z)`
    (d) `f(z)`

## Case study: numerical integration

To conclude this chapter, we will develop a simple numerical integration tool, and along the way, illustrate the use of many properties of first-class functions. Each step is driven by a desire to make our approach more general and to reduce duplication.  The idea behind numerical integration is simple: we want to find the area under the curve by approximating a complex curve with simpler components.

The two simpliest approaches are the __midpoint__ and __trapezoid__ rules; the mid point rule approximates a curve by a rectangle, and the trapezoid rule by a trapezoid. Each takes a function we want to integrate, `f`, and a range to integrate over, from `a` to `b`. For this example we'll try to integrate `sin x` from 0 to pi, because it has a simple answer: 2.

```R
midpoint <- function(f, a, b) {
  (b - a) * f((a + b) / 2)
}

trapezoid <- function(f, a, b) {
  (b - a) / 2 * (f(a) + f(b))
}

midpoint(sin, 0, pi)
trapezoid(sin, 0, pi)
```

Neither of these functions gives a very good approximation, so we'll do what we normally do in calculus: break up the range into smaller pieces and integrate each piece using one of the simple rules. This is called __composite integration__, and we'll implement it with two new functions:

```R
midpoint_composite <- function(f, a, b, n = 10) {
  points <- seq(a, b, length = n + 1)
  h <- (b - a) / n
  
  area <- 0
  for (i in seq_len(n)) {
    area <- area + h * f((points[i] + points[i + 1]) / 2)
  }
  area
}

trapezoid_composite <- function(f, a, b, n = 10) {
  points <- seq(a, b, length = n + 1)
  h <- (b - a) / n
  
  area <- 0
  for (i in seq_len(n)) {
    area <- area + h / 2 * (f(points[i]) + f(points[i + 1]))
  }
  area
}

midpoint_composite(sin, 0, pi, n = 10)
midpoint_composite(sin, 0, pi, n = 100)
trapezoid_composite(sin, 0, pi, n = 10)
trapezoid_composite(sin, 0, pi, n = 100)
    
mid <- sapply(1:20, function(n) midpoint_composite(sin, 0, pi, n))
trap <- sapply(1:20, function(n) trapezoid_composite(sin, 0, pi, n))
matplot(cbind(mid = mid, trap))
```

But notice that there's a lot of duplication across `midpoint_composite` and `trapezoid_composite`: they are basically the same apart from the internal rule used to integrate over a simple range. Let's extract out a general composite integrate function:

```R
composite <- function(f, a, b, n = 10, rule) {
  points <- seq(a, b, length = n + 1)
  
  area <- 0
  for (i in seq_len(n)) {
    area <- area + rule(f, points[i], points[i + 1])
  }
  
  area
}

midpoint_composite(sin, 0, pi, n = 10)
composite(sin, 0, pi, n = 10, rule = midpoint)
composite(sin, 0, pi, n = 10, rule = trapezoid)
```

This function now takes two functions as arguments: the function to integrate, and the integration rule to use for simple ranges. We can now add even better rules for integrating small ranges:

```R
simpson <- function(f, a, b) {
  (b - a) / 6 * (f(a) + 4 * f((a + b) / 2) + f(b))
}

boole <- function(f, a, b) {
  pos <- function(i) a + i * (b - a) / 4
  fi <- function(i) f(pos(i))
  
  (b - a) / 90 * 
    (7 * fi(0) + 32 * fi(1) + 12 * fi(2) + 32 * fi(3) + 7 * fi(4))
}
```
    
Let's compare these different approaches.

```R
expt1 <- expand.grid(
  n = 5:50, 
  rule = c("midpoint", "trapezoid", "simpson", "boole"), 
  stringsAsFactors = F)

abs_sin <- function(x) abs(sin(x))
run_expt <- function(n, rule) {
  composite(abs_sin, 0, 4 * pi, n = n, rule = match.fun(rule))
}

library(plyr)
res1 <- mdply(expt1, run_expt)

library(ggplot2)
qplot(n, V1, data = res1, colour = rule, geom = "line")
```

It turns out that the midpoint, trapezoid, Simpson and Boole rules are all examples of a more general family called Newton-Cotes rules. (They are polynomials of increasing complexity).  We can take our integration one step further by extracting out this commonality to produce a function that can generate any general Newton-Cotes rule:

```R
# http://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas
newton_cotes <- function(coef, open = FALSE) {
  n <- length(coef) + open
  
  function(f, a, b) {
    pos <- function(i) a + i * (b - a) / n
    points <- pos(seq.int(0, length(coef) - 1))
    
    (b - a) / sum(coef) * sum(f(points) * coef)        
  }
}

trapezoid <- newton_cotes(c(1, 1))
midpoint <- newton_cotes(1, open = TRUE)
simpson <- newton_cotes(c(1, 4, 1))
boole <- newton_cotes(c(7, 32, 12, 32, 7))
milne <- newton_cotes(c(2, -1, 2), open = TRUE)

expt1 <- expand.grid(n = 5:50, rule = names(rules), 
  stringsAsFactors = FALSE)
run_expt <- function(n, rule) {
  composite(abs_sin, 0, 4 * pi, n = n, rule = rules[[rule]])
}
```

Mathematically, the next step in improving numerical integration is to move from a grid of evenly spaced points to a grid where the points are closer together near the end of the range, such as __Gaussian quadrature__.  That's beyond the scope of this case study, but you would use similar techniques to add it.

### Exercises

* Instead of creating individual fuctions `midpoint()`, `trapezoid()`, `simpson()` etc, we could store them in a list. If we do that, how does the code change? Can you create the list of functions from a list of coefficients for the Newton-Cotes formulae?

* The tradeoff in integration rules is that more complex rules are slower to compute, but need fewer pieces.  For `sin()` in the range [0, pi], determine the number of pieces needed to for each rule to be equally accurate. Illustrate your results with a graph.  How do they change for different functions? `sin(1 / x^2)` is particularly challenging. 

*  For each of the Newton-Cotes rules, how many pieces do you need to get within 0.1% of the true answer for `sin()` in the range [0, pi].  Write a function that determines that automatically for any function (hint: look at `optim()` and construct a one-argument function with closures)