```{r setup, include=FALSE}
knitr::opts_chunk$set(comment = "", echo = TRUE)
```

## Packages to run this presentation

```{r packages, message=FALSE, warning=FALSE}
library("readr")
library("tibble")
library("dplyr")
library("tidyr")
library("stringr")
library("ggplot2")
library("purrr")
library("broom")
```

## Motivation

As you know, purrr is a recent package from Hadley Wickham, focused on lists and functional programming, like dplyr is focused on data-frames. 

I figure a good way to learn a new package is to try to solve a problem, so we have a dataset: 

- you can [view](https://github.com/ijlyttle/isugg_purrr/blob/gh-pages/temperature.csv) or [download](http://ijlyttle.github.io/isugg_purrr/temperature.csv)

- you can download the [source](http://ijlyttle.github.io/isugg_purrr/presentation.Rmd) of this presentation

- these are three temperatures recorded simultaneously in a piece of electronics

- it will be very valuable to be able to characterize the transient temperature for each sensor

- we want to apply the same set of models across all three sensors

- it will be easier to show using pictures

## Let's get the data into shape

Using the readr package

```{r load}
temperature_wide <- 
  read_csv("temperature.csv") %>%
  print()
```

## Is `temperature_wide` "tidy"?

```{r is_tidy, echo=FALSE}
print(temperature_wide)
```

Why or why not?

## Tidy data

1. Each column is a variable
2. Each row is an observation
3. Each cell is a value

(http://www.jstatsoft.org/v59/i10/paper)

My personal observation is that "tidy" can depend on the context, on what you want to do with the data.

## Let's get this into a tidy form

```{r}
temperature_tall <-
  temperature_wide %>%
  gather(key = "id_sensor", value = "temperature", starts_with("temp")) %>%
  mutate(id_sensor = str_replace(id_sensor, "temperature_", "")) %>%
  print()
```

## Now, it's easier to visualize

```{r}
temperature_tall %>%
  ggplot(aes(x = instant, y = temperature, color = id_sensor)) +
  geom_line()
```

## Rearrange a bit more

**`delta_time`** $\Delta t$ 

change in time since event started, s

**`delta_temperature`**: $\Delta T$

change in temperature since event started, °C

```{r}
delta <- 
  temperature_tall %>%
  arrange(id_sensor, instant) %>%
  group_by(id_sensor) %>%
  mutate(
    delta_time = as.numeric(instant) - as.numeric(instant[[1]]),
    delta_temperature = temperature - temperature[[1]]
  ) %>%
  select(id_sensor, delta_time, delta_temperature)
```

## Let's have a look

```{r}
delta %>%
  ggplot(aes(x = delta_time, y = delta_temperature, color = id_sensor)) +
  geom_line()  
```

## Curve-fitting

We want to see how three different curve-fits might perform on these three data-sets:

### Newtonian cooling

\[\Delta T = \Delta {T_0} \left[ 1 - \exp \left( { - \frac}} \right) \right] \]

### Semi-infinite solid

\[\Delta T = \Delta {T_0}\operatorname{erfc} \left( {\sqrt {\frac}} } \right)\]

### Semi-infinite solid with convection

\[\Delta T = \Delta {T_0}\left[ {\operatorname{erfc} \left( {\sqrt {\frac}} } \right) - \exp \left( {B{i_0} + \frac{4}\frac}} \right)\operatorname{erfc} \left( {\sqrt {\frac}}  + \frac}{2}\sqrt {\frac}} } \right)} \right]\]

## Some definitions

```{r}
# reference: http://stackoverflow.com/questions/29067916/r-error-function-erfz
# (see Abramowitz and Stegun 29.2.29)
erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1
erfc <- function(x) 2 * pnorm(x * sqrt(2), lower = FALSE)
```

```{r}
newton_cooling <- function(x) {
  nls(
    delta_temperature ~ delta_temperature_0*(1 - exp(-delta_time/tau_0)),
    start = list(delta_temperature_0 = -10, tau_0 = 50),
    data = x
  )
}
```

## More math

```{r}
semi_infinite_simple <- function(x) {
  nls(
    delta_temperature ~ delta_temperature_0*erfc(sqrt(tau_0/delta_time)),
    start = list(delta_temperature_0 = -10, tau_0 = 50),
    data = x
  )    
}

semi_infinite_convection <- function(x){
  nls(
    delta_temperature ~
      delta_temperature_0*(
        erfc(sqrt(tau_0/delta_time)) -
        exp(Bi_0 + (Bi_0/2)^2*delta_time/tau_0)*
          erfc(sqrt(tau_0/delta_time) + 
        (Bi_0/2)*sqrt(delta_time/tau_0))
      ),
    start = list(delta_temperature_0 = -5, tau_0 = 50, Bi_0 = 1.e6),
    data = x
  )
}
```

## Before we get into purrr

Before doing anything, we want to show that we can do something with one dataset and one model-function:

```{r}
tmp_data <- delta %>% filter(id_sensor == "a")

tmp_model <- newton_cooling(tmp_data)

summary(tmp_model)
```

## Look at predictions

```{r}
tmp_pred <- 
  tmp_data %>%
  mutate(modeled = predict(tmp_model, data = .)) %>%
  select(id_sensor, delta_time, measured = delta_temperature, modeled) %>%
  gather("type", "delta_temperature", measured:modeled) %>%
  print()
```

## A more-useful look

```{r}
tmp_pred %>%
  ggplot(aes(x = delta_time, y = delta_temperature, linetype = type)) +
  geom_line()
```

## "Regular" data-frame

```{r}
print(delta)
```

Each column of the dataframe is a vector - in this case, a character vector and two doubles

## How to make a weird data-frame

Here's where the fun starts - a column of a data-frame can be a list.

- use `tidyr::nest()` to makes a column `data`, which is a list of data-frames

- this seems like a stronger expression of the `dplyr::group_by()` idea 

```{r}
delta_nested <- 
  delta %>%
  nest(-id_sensor) %>%
  print()
```

## Map data-frames to the modeling function

- `map()` is like `lapply()`

- `map()` returns a list-column (it keeps the weirdness)

```{r}
model_nested <-
  delta_nested %>%
  mutate(model = map(data, newton_cooling)) %>%
  print()
```

## We can use `map2()` to make the predictions

- `map2()` is like `mapply()`

- designed to map two colunms (`model`, `data`) to a function `predict()`

```{r}
predict_nested <-
  model_nested %>%
  mutate(pred = map2(model, data, predict)) %>%
  print()
```

## We need to get out of the weirdness

- use `unnest()` to get back to a regular data-frame

```{r}
predict_unnested <- 
  predict_nested %>%
  unnest(data, pred) %>% 
  print()
```

## We can wrangle the predictions

- get into a form that makes it easier to plot

```{r}
predict_tall <- 
  predict_unnested %>%
  rename(modeled = pred, measured = delta_temperature) %>%
  gather("type", "delta_temperature", modeled, measured) %>%
  print()
```

##  We can visualize the predictions

```{r}
predict_tall %>%
  ggplot(aes(x = delta_time, y = delta_temperature)) +
  geom_line(aes(color = id_sensor, linetype = type))
```

## Now we want to look at a selection of models

Make a list of functions to model:

```{r}
list_model <-
  list(
    newton_cooling = newton_cooling,
    semi_infinite_simple = semi_infinite_simple,
    semi_infinite_convection = semi_infinite_convection
  )
```

## Step: write a function to define the "inner" loop 

```{r}
fn_model <- function(.model, df){
  # safer to avoid non-standard evaluation
  # df %>% mutate(model = map(data, .model)) 
  
  df$model <- map(df$data, possibly(.model, NULL))
  df
}
```

- for a given model-function and a given (weird) data-frame, return a modified version of that data-frame with a column `model`, which is the model-function applied to each element of the data-frame's `data` column (which is itself a list of data-frames)

- the purrr functions `safely()` and `possibly()` are **very** interesting. I think they could be useful outside of purrr as a friendlier way to do error-handling.

## Step: `map_df()` to define the "outer" loop

```{r}
model_nested_new <-
  list_model %>%
  map_df(fn_model, delta_nested, .id = "id_model") %>%
  print()
```

- for each element of a list of model-functions, run the inner-loop function, and row-bind the results into a data-frame

- we want to discard the rows where the model failed
- we also want to investigate why they failed, but that's a different talk

## Step: `map()` to identify the null models

```{r}
model_nested_new <-
  list_model %>%
  map_df(fn_model, delta_nested, .id = "id_model") %>%
  mutate(is_null = map(model, is.null)) %>%
  print()
```

- using `map(model, is.null)` returns a list column
- to use `filter()`, we have to escape the weirdness

## Step: `map_lgl()` to identify nulls and get out of the weirdness

```{r}
model_nested_new <-
  list_model %>%
  map_df(fn_model, delta_nested, .id = "id_model") %>%
  mutate(is_null = map_lgl(model, is.null)) %>%
  print()
```

- using `map_lgl(model, is.null)` returns a vector column

## Step: `filter()` and `select()` to clean up

```{r}
model_nested_new <-
  list_model %>%
  map_df(fn_model, delta_nested, .id = "id_model") %>%
  mutate(is_null = map_lgl(model, is.null)) %>%
  filter(!is_null) %>%
  select(-is_null) %>%
  print()
```

## Let's get predictions

```{r}
predict_nested <- 
  model_nested_new %>%
  mutate(pred = map2(model, data, predict)) %>%
  print()
```

## `unnest()`, make it tall

```{r}
predict_tall <-
  predict_nested %>%
  unnest(data, pred) %>% 
  rename(modeled = pred, measured = delta_temperature) %>%
  gather("type", "delta_temperature", modeled, measured) %>%
  print()
```

##  We can visualize the predictions

```{r}
predict_tall %>%
  ggplot(aes(x = delta_time, y = delta_temperature)) +
  geom_line(aes(color = id_sensor, linetype = type)) +
  facet_grid(id_model ~ .)
```

## Let's get the residuals

```{r}
resid <-
  model_nested_new %>%
  mutate(resid = map(model, resid)) %>%
  unnest(data, resid) %>%
  print()
```

## And visualize them

```{r}
resid %>%
  ggplot(aes(x = delta_time, y = resid)) +
  geom_line(aes(color = id_sensor)) +
  facet_grid(id_model ~ .)
```

## Using broom package to look at model-statistics

The `tidy()` function extracts statistics from a model

```{r}
model_parameters <- 
  model_nested_new %>%
  select(id_model, id_sensor, model) %>%
  mutate(tidy = map(model, tidy)) %>%
  select(-model) %>%
  unnest() %>%
  print()
```

## Get a sense of the coefficients

```{r}
model_summary <-
  model_parameters %>%
  select(id_model, id_sensor, term, estimate) %>%
  spread(key = "term", value = "estimate") %>%
  print()
```

## Summary

- this is just a smalll part of purrr
- there seem to be parallels between `tidyr::nest()/purrr::map()` and `dplyr::group_by()/dplyr::do()`
    - to my mind, the purrr framework is more understandable
    - update tweet from [Hadley](https://twitter.com/hadleywickham/status/719542847045636096)
    
References from Hadley:

- [purrr 0.1.0 announcement](http://blog.rstudio.org/2015/09/29/purrr-0-1-0/)
- [purrr 0.2.0 announcement](http://blog.rstudio.org/2016/01/06/purrr-0-2-0/)
- [chapter from Garrett Grolemund and Hadley's forthcoming book](http://r4ds.had.co.nz/iteration.html)  
