---
title: "Minimal Advice on Programming in R"
author: 36-402, Section A
date: 29 January 2019 (Lecture 5)
output: slidy_presentation
---

```{r}
library(knitr)
opts_chunk$set(size="small", cache=TRUE,
               autodep=TRUE, warnings=TRUE, error=TRUE)
```

## Why do _you_ need to do programming?

> - Programs for data analysis are tools or machines
>    + Understand your tools!
> - If you know how to make new tools, you always have the right tools
>    + Don't distort your problem to match old tools
> - Programming is making those tools
> - Programming is also making your ideas clear

## Agenda for today

1. Style
2. Functions
3. Debugging
4. Avoiding iteration
5. Design advice


## Style: Why?

> - Code is about telling the computer what to do
> - Code is _also_ about communicating with other people
>     + Including communicating with _yourself_, later


## Style: Commandments

- Use meaningful names for variables, data-sets, functions, etc.
    + Avoid magic numbers
- Indent your code blocks
- **Comment your code**
	+ Comment everything that seems clever or tricky
    + Comment every trick you took from somewhere else
    + Comment every function with purpose, inputs, outputs


## Functions

- Use functions to tie together related commands so they can be repeated
    + Less work (in the long run)
	+ Fewer errors (get it right _once_, make improvements _once_)
	+ Clearer: allows for abstraction
	+ Lets you build complicated procedures out of modular, re-usable parts


## Functions: Example with Pareto distribution

- Distribution of wealth has very, very heavy tails
- A good (but not perfect) model is the **Pareto** distribution, pdf
\[
f(x;\alpha,x_0) = \begin{cases}
\frac{\alpha-1}{x_0}\left(\frac{x}{x_0}\right)^{-\alpha} & x \ge x_0\\
0 & x < x_0
\end{cases}
\]
- We often want the quantiles (where does the 99% begin? the 99.9%? the 80%?)
    + The **mathematical** quantile function is
\[
Q(p; \alpha,x_0) = x_0(1-p)^{-\frac{1}{\alpha-1}}
\]
- Could do stuff like
```{r}
# alpha = 2.33, x_0 = 6e8, p = 0.5
6e8 * (1-0.5)^(-1/(2.33-1))

# alpha = 2.33, x_0 = 6e8, p = 0.4
6e8 * (1-0.4)^(-1/(2.33-1))

# alpha = 2.5, x_0 = 1e6, p = 0.92
1e6 * (1-0.92)^(-1/(2.5-1))
```

## Functions: Pareto example

- Bundle up the repetitive code into one function:

```{r}
# Calculate quantiles of the Pareto distribution
# Inputs: desired quantile (p)
  # exponent of the distribution (exponent)
  # lower threshold of the distribution (threshold)
# Outputs: the pth quantile
qpareto.1 <- function(p, exponent, threshold) {
  q <- threshold*((1-p)^(-1/(exponent-1)))
  return(q)
}
```
- Now we can say:
```{r}
qpareto.1(p=0.5, exponent=2.33, threshold=6e8)
qpareto.1(p=0.4, exponent=2.33, threshold=6e8)
qpareto.1(p=0.92, exponent=2.5, threshold=1e6)
```



## Parts of a function

- Arguments: what does the function need to know?
    + Default values let you omit arguments that don't change often
    + Default values can always be over-ridden
	+ [No-argument functions are possible](#noargs)
- Return value: what does the function give back to R?
    + Try to be explicit about what you return
    + Use `c()` or `list()` to return more than one thing (We will see some examples of this)
- The body: do these things to the arguments to produce the return value

## Parts of a function: Pareto quantiles cont'd

```{r}
# Calculate quantiles of the Pareto distribution
# Inputs: desired quantile (p)
  # exponent of the distribution (exponent)
  # lower threshold of the distribution (threshold)
  # flag for whether to give lower or upper quantiles (lower.tail)
# Outputs: the pth quantile
qpareto.2 <- function(p, exponent, threshold, lower.tail=TRUE) {
  if(lower.tail==FALSE) {
    p <- 1-p
  }
  q <- threshold*((1-p)^(-1/(exponent-1)))
  return(q)
}
```

```{r}
qpareto.2(p=0.92, exponent=2.5, threshold=1e6)
qpareto.2(p=0.92, exponent=2.5, threshold=1e6, lower.tail=FALSE)
```

## When should we write functions?

> - When we're repeating the same thing multiple times
>     + Especially when: we're repeating **almost** the same thing
> - When we're writing very similar code over and over
> - When we're breaking down a big job into smaller ones
> - When we're told to

## Chaining functions together

- Divide the big job ("come up with a confidence band") into smaller parts
- Write a function for each part
- Make sure the outputs of the early steps can be used as inputs to the later steps
- Repeat (recurse) if need be: the sub-functions spawn sub-sub-functions, etc.
- We can often re-use the sub-functions

## Calling functions from other functions

- Functions can call other functions, including other functions we've written
- Could have avoided replicating the quantile code from `qpareto.1` in `qpareto.2`:
```{r}
# Calculate quantiles of the Pareto distribution
# Inputs: desired quantile (p)
  # exponent of the distribution (exponent)
  # lower threshold of the distribution (threshold)
  # flag for whether to give lower or upper quantiles (lower.tail)
# Outputs: the pth quantile
qpareto.3 <- function(p, exponent, threshold, lower.tail=TRUE) {
  if(lower.tail==FALSE) {
    p <- 1-p
  }
  # Call qpareto.1() to do the actual calculation
  q <- qpareto.1(p, exponent, threshold)
  return(q)
}
```


## Functions calling other functions (cont'd)

- You will often be working with several intertwined functions at the same time
- The called functions have to be defined when they are _used_
    + Example: the lecture code for `cv.lm()` calls a function, `response.name()`,
	to figure out what variable in a `lm` formula is the response
	+ `response.name()` is defined just after `cv.lm()`, but before the code +runs_ `cv.lm()`
- If you change (fix) the definition of a function, the new, revised one will
be used thereafter (no snapshots)



## Debugging

> - Eventually, you will break something.
> - **Problem**: the final piece of code you write gives an incomprehensible error message
> - **Idea**: Carefully characterize the bug, then narrow down where it could be until the the problem is obvious.
> - Strategy: make sure all the lower-level pieces are working properly
>     + Start debugging from the bottom
>     + If all the small bits work on their own, the bug comes when you put them together

## Debugging: Test cases

> - Use test cases
>     + "What should this do if it works?"
>     + "What should this _not_ do if it's broken?"
> - When you meet a new bug, write a new test
>     + "How do I know if this is fixed?"
>     + Test could be very minimal: "It runs without crashing"

## Debugging: Some tools

- `traceback()`: Show all the calls leading to the last error
    + Look for the last line of code you wrote
    + This helps localize where things went wrong
- `print()`: Make R send output to the console
    + Most useful a bit before where things crash/go crazy
	+ Remember to get rid of it when the bug is fixed!
- `browser()`: A convenient way to mess around inside a function
    + When called within a function, gives you an R shell there
    + Avoids problems with local variables disappearing
    + A common trick: Executes `[Dangerous Expression]` and calls `browser()` if there is an error.
```
tryCatch({[Dangerous Expression]},warning=function(w){browser()},
         error=function(e){browser()})
```

## Debugging: Pareto example continued

- I want to generate a  Pareto random variable
    + If $F(x)$ is the CDF for a distribution and $U$ is a uniform random variable, then $F^{-1}(U)$ has the $F$ distribution.

- New version of the quantile function that sanity-checks the arguments:
```{r}
# Calculate quantiles of the Pareto distribution
# Inputs: desired quantile (p)
  # exponent of the distribution (exponent)
  # lower threshold of the distribution (threshold)
  # flag for whether to give lower or upper quantiles (lower.tail)
# Outputs: the pth quantile
qpareto.4 <- function(p, exponent, threshold, lower.tail=TRUE) {
  stopifnot(p >= 0, p <= 1, exponent > 1, threshold > 0)
  q <- qpareto.3(p,exponent,threshold,lower.tail)
  return(q)
}
```

- Try it with a bad input
```{r}
qpareto.4(p=0.92,exponent=2.5,threshold=-1,lower.tail=FALSE)
```

## Debugging: Pareto example continued

- Now define a random-Pareto function.  This version is (intentionally) buggy!
```{r}
# Generate random numbers from the Pareto distribution
# Inputs: number of random draws (n)
  # exponent of the distribution (exponent)
  # lower threshold of the distribution (threshold)
# Outputs: vector of random numbers
rpareto <- function(n,exponent,threshold) {
  x <- vector(length=n)
  for (i in 1:n) {
    x[i] <- qpareto.4(p=rnorm(1),exponent=exponent,threshold=threshold)
  }
  return(x)
}
```

## Debugging: Pareto example continued
- Our buggy code gives us errors!
```{r}
rpareto(10)
```

```{r}
rpareto(n=10,exponent=2.5,threshold=1)
```

## Debugging: The in-class exercise {.smaller}

EXERCISE: Where is the bug?  Why doesn't `rpareto(n=10,exponent=2.5,threshold=1)` work?

```{r, eval=FALSE}
qpareto.1 <- function(p, exponent, threshold) {
  q <- threshold*((1-p)^(-1/(exponent-1)))
  return(q)
}

qpareto.3 <- function(p, exponent, threshold, lower.tail=TRUE) {
  if(lower.tail==FALSE) {
    p <- 1-p
  }
  q <- qpareto.1(p, exponent, threshold)
  return(q)
}

qpareto.4 <- function(p, exponent, threshold, lower.tail=TRUE) {
  stopifnot(p >= 0, p <= 1, exponent > 1, threshold > 0)
  q <- qpareto.3(p,exponent,threshold,lower.tail)
  return(q)
}

rpareto <- function(n,exponent,threshold) {
  x <- vector(length=n)
  for (i in 1:n) {
    x[i] <- qpareto.4(p=rnorm(1),exponent=exponent,threshold=threshold)
  }
  return(x)
}
```

## Exercise: Solution

- We've checked the `qpareto` functions
- The error message is that `p >= 0` is not true
- The error message comes from `qpareto.4()`
- So a negative value of `p` has to be going into `qpareto.4()` when it's called
by `rpareto()`
- And in fact that code says `p=rnorm(1)`, and `rnorm()` draws from the standard normal, which is negative half the time
- That should be `p=runif(1)` (uniform distribution on $[0,1]$), not `p=rnorm(1)`
    + As mentioned, this is the bug I actually made when I first wrote `rpareto` in 2004


## Avoid iteration

> - Iterating in R is often quite slow
>    + Not (always) true in other languages
>    + Basically: every time you assign in R, you make a copy
> - R has lots of ways to avoid iteration
> - This is often faster!  And easier to follow!

## Avoid iteration: Pareto example cont'd

- A very iterative way of making many random Pareto draws:
```{r}
# Generate random numbers from the Pareto distribution
# Inputs: number of random draws (n)
  # exponent of the distribution (exponent)
  # lower threshold of the distribution (threshold)
# Outputs: vector of random numbers
rpareto.1 <- function(n,exponent,threshold) {
  x <- c()
  for (i in 1:n) {
    x <- c(x,qpareto.4(p=runif(1),exponent=exponent,threshold=threshold))
  }
  return(x)
}
```

```{r}
system.time(rpareto.1(n=1e4,exponent=2.5,threshold=1))
```

## Avoid iteration: Pareto example cont'd

- A not-so-iterative way:
```{r}
# Generate random numbers from the Pareto distribution
# Inputs: number of random draws (n)
  # exponent of the distribution (exponent)
  # lower threshold of the distribution (threshold)
# Outputs: vector of random numbers
rpareto.2 <- function(n,exponent,threshold) {
  x = replicate(n, qpareto.4(p=runif(1),exponent=exponent,threshold=threshold))
  return(x)
}


system.time(rpareto.2(n=1e4,exponent=2.5,threshold=1))
```

## Avoiding iteration: Vectorizing

- Lots of R functions are happy with _vector_ arguments, and **vectorize** along them

- Our `qpareto()` functions can handle vector `p` arguments just fine:
```{r}
qpareto.4(c(.1,.2,.3), 2.5, 1)
```

- So use that:
```{r}
# Generate random numbers from the Pareto distribution
# Inputs: number of random draws (n)
  # exponent of the distribution (exponent)
  # lower threshold of the distribution (threshold)
# Outputs: vector of random numbers
rpareto.3 <- function(n,exponent,threshold) {
  x = qpareto.4(p=runif(n), exponent=exponent, threshold=threshold)
  return(x)
}
```

```{r}
system.time(rpareto.3(10000,exponent=2.5,threshold=1))
```

## Avoiding iteration: `ifelse`

- use `ifelse` to do conditional things to entire vectors

```{r, fig.width=4, fig.height=4}
# Calculate Huber's loss function
# Input: vector of numbers x
# Return: x^2 for |x|<1, 2|x|-1 otherwise
huber <- function(x) {
  return(ifelse(abs(x) <= 1, x^2, 2*abs(x)-1))
}
x = seq(from=-3,to=3,length.out=100)
plot(x,huber(x), type="l")
```


## Avoiding iteration: matrix and vector operations

```{r}
# Get some long vectors
n = 1e6
a = rnorm(n)
b = rnorm(n)
# Find the inner product with a for() loop
system.time({
  total = 0
  for(i in 1:length(a)){total = total + a[i]*b[i]}
})
# Find the inner product with matrix operations
system.time({
a %*% b
})
```


## Avoiding iteration: `apply` and friends

- Use `apply` to apply a function to rows or columns of a matrix
- Use `sapply` to apply a function to every element of a vector
```{r}
#Make a 100x100 matrix of uniform [0,1] variables
X = matrix(runif(10000),nrow=100,ncol=100)
#Find the maximum of each column
column.maxs = apply(X, MARGIN = 2, FUN = max)
#Plot a histogram.  Do you know the distribution?
hist(column.maxs,20)
```

- The `tidyverse` package defines a _lot_ of new functions to push this idea to logical extremes


## Do not be afraid of iteration

- Some jobs are just inherently iterative
- Use it when that's what you need to do
- If you're storing stuff as you go, create the variables which you're storing into
outside the loop
```{r, eval=FALSE}
sim.ar <- function(n,x1,phi=0,sd=1) {
  x <- vector(length=n)
  x[1] <- x1
  for (t in 2:n) {
    x[t] <- phi*x[t-1]+rnorm(1,mean=0,sd=sd)
  }
  return(x)
}
```

## Returning vectors and lists

- A function can return one *thing*.  Often you want to return more than one value:
    + For several values of the same type, you can return a vector
    +  For several different sorts of values, you can return a list.  You can even assign names!
- Returning a vector: Look at `rpareto()`:
```{r, eval=FALSE}
rpareto.3 <- function(n,exponent,threshold) {
  x = qpareto.4(p=runif(n), exponent=exponent, threshold=threshold)
  return(x) #If n>1, this is a vector!
}
```
- Returning a list: look at the end of `cv_bws_npreg()` from lecture 4:
```{r, eval=FALSE}
  return(list(best.bw=best.bw,CV_MSEs=CV_MSEs,
              fold_MSEs=fold_MSEs))
```

## Design advice

- Design from the top down
    + Start from _what_ you want to do (not how)
	+ Break that task into a few (2--6) smaller tasks
	+ Keep breaking tasks into smaller tasks until you can't go further
- Code from the bottom up
    + Start with the smallest, simplest tasks
	+ Write a function for each one, and test it
	+ Then write functions that call _those_ functions (and test them)
	+ Keep going until you have top-level functions for the whole task

## What you need to remember

1. **Use meaningful names** and **Comment your code**
2. Write functions to bundle related commands together
3. Split big jobs into a bunch of related functions
4. Practice debugging
5. Avoid iteration, in favor of acting on whole objects


## Aside: Style: Things not to do with your main data frame

- Don't use `attach()`
    + Makes it really hard to track down errors
	+ Can shorten your code using `with()` instead
- Try not to call your main data set `data`
    - vague
    - the name of a function for loading data sets from packages
    - often the name of arguments in  model-fitting functions
    - both of those make it hard to trace what's going on
```{r,eval=FALSE}
# Not what you want to read:
data <- read.csv("http://www.stat.cmu.edu/~cshalizi/uADA/19/hw/02/uval.csv")
data.lm <- lm(growth~log(gdp)+underval, data=data)
```




## Aside: No-argument functions {#noargs}

- No-argument functions can be useful for simulation

```{r, fig.height=3,fig.width=4}
# Function to return the biggest correlation of 5 length-10 Gaussian vectors
# Takes no arguments, so it can be used in a simulation
random.correlation.max = function(){
  #Draw a 10x5 matrix of Gaussians
  X = matrix(rnorm(50),nrow=10,ncol=5)
  #Compute the correlation matrix
  cor.matrix = cor(X)
  #Find the biggest entry, making sure to avoid the diagonal (which is all 1)
  return(max(cor.matrix[upper.tri(cor.matrix)]))
}
```
- Draw 100 samples and plot them

```{r}
hist(replicate(100,random.correlation.max()),20,col='grey', main="Max. Correlation",xlab='correlation')
```


## Aside: Functions refer to the _current_ definition of other functions

- A very simple function that adds one to its argument

```{r}
add.one = function(x){
  x+1
}
```
- A function that calls `add.one()` to do its addition
```{r}
add.stuff = function(x){
  return(add.one(x))
}
```
```{r}
add.stuff(3)
```
- A terrible job rewriting `add.one()`:
```{r}
add.one = function(x){
  x+2
}
```
- Now what?
```{r}
add.stuff(3)
```

