---
title: "The PEcAn RTM package"
author: Alexey Shiklomanov
output_format: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{The PEcAn RTM package}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r, setup, include=FALSE}
knitr::opts_chunk$set(fig.width = 6, fig.height = 4.5, out.width = '100%', dpi=60)
```

# Introduction

The PEcAnRTM package provides tools for analyses involving common radiative transfer models. The highlights of this package are its ability to efficiently run a suite of related leaf and canopy radiative transfer models, as well as to perform maximum likelihood and, particularly, Bayesian inversions of these models. 

# Installation
The easiest way to install this package is via `install_github` from the `devtools` package.

```{r, eval=FALSE}
install.packages("devtools")
library(devtools)
install_github("ashiklom/pecan", subdir = "modules/rtm")
# Defaults to branch 'master'. 
# For custom branches, add `ref = "branchname"`
```

This package relies on a modern (>= 2003) Fortran compiler, as determined by your R installation. On most Unix systems, this is standard, but R may specify a particular compiler version that you don't have, resulting in an installation error. To fix this, simply add the following to your system `~/R/Makevars` file.

```
FC = gfortran
```

# Overview of features

## Simulating reflectance

Available radiative transfer models are called by passing a vector of their parameters. Similar models with different versions (e.g. PROSPECT 4, 5, 5B) also have a "version" argument. These models return a matrix of reflectance, transmittance, and/or absorption spectra (depending on the model) at 1 nm resolution over the wavelength range 400 to 2500 nm.

The PROSPECT family of models returns the reflectance (column 1) and transmittance (column 2) for an individual leaf as a function of 4 to 6 parameters (depending on the version).

```{r}
library(PEcAnRTM)
wl <- 400:2500
params <- c(
  "N" = 1.4, "Cab" = 40, "Car" = 15,
  "Cbrown" = 0.5, "Cw" = 0.002, "Cm" = 0.004
)

p4 <- prospect(params[c(-3, -4)], version = 4)
p4
```

Notice only the first and last few lines of the spectra are printed, and the wavelength is annotated on the side.
This is is because all RTM simulations in this package are of a special matrix class `spectra`.

```{r}
class(p4)
```

This class allows traditional subsetting via index with single brackets (`[]`), as well as subsetting by wavelength with double brackets.

```{r}
p4[1:50, 1]
p4[[500:520, 2]]
```

The package also provides special plotting (`plot`, `matplot`) and combining (`cbind`) methods. Note that the `cbind` method automatically matches wavelengths, which facilitates working with spectra with different wavelengths.

```{r}
p5 <- prospect(params[-4], version = 5)
p5b <- prospect(params, version = "5B")
p_multi <- cbind(p4, p5, p5b)
matplot(p_multi[, c(1, 3, 5)], lty = 1:3, col = 2, ylim = c(0, 1))
matplot(1 - p_multi[, c(2, 4, 6)], lty = 1:3, col = 3,  add = TRUE)
legend("topright", c("Reflectance", "Transmittance"), col=c(2, 3), lty = 1)
legend("top", c("4", "5", "5B"), lty = 1:3)
```

The SAIL family of models returns the bidirectional (1), hemispherical directional (2), directional hemispherical (3), and bidirectional hemispherical (4) reflectance factors for a canopy with a given set of approximately 20 parameters. It is often coupled to the PROSPECT model as PRO4SAIL. (Again, note that the return type is a `spectra`, which leads `matplot` to automatically put wavelengths on the x axis.)

```{r}
sail.params <- defparam("pro4sail")
print(sail.params)
p4s <- pro4sail(sail.params)
matplot(p4s, xlab="Wavelength (nm)", ylab="Reflectance")
legend("topright", as.character(1:4), col=1:4, lty=1:4)
```

The above example illustrates the use of `defparam` to get the default parameters for a particular model. Similarly, `model.list` is a `data.frame` containing all currently available models.

```{r}
print(model.list)
```


## Inversion

A novel feature of this package is the ability to perform a Bayesian inversion of a Radiative Transfer Model. Here are several advantages of the Bayesian approach:

* **Parameter uncertainty:** The output of a Bayesian analysis is a full joint probability distribution of the model parameters, which includes a robust estimate of their uncertainy and covariance between parameters.
* **Prior knowledge:** If previous, independent estimates of parameters are available, these can be used to inform the model.
* **Partitioning variability:** Random effects models provide a powerful framework for understanding the sources of variability and uncertainty in a data set.

An inversion can be performed using the `invert.auto` function, which uses the Metropolis Hastings MCMC algorithm to invert an arbitrary model. There are a lot of configuration options to `invert.auto`, so the recommended way to perform an inversion is to start with a default settings list, provided by the package itself. In the following sample, we demonstrate the default inversion of the PROSPECT model.

```{r}
invert.options <- default.settings.prospect
```

The model to be inverted is, in this case, just a one-line call to the PROSPECT 5 model with the params vector. It returns a vector of reflectance values. The requirement is that the "model" output be as long as the length (or number of rows) of the observation vector (or matrix). In this case, the PROSPECT model returns 2101 reflectance values, so our observation also has to have this many points.

```{r}
invert.options$model
```

The recommended way to set settings is to start with a default object and modify it. For example, in the following block, we reduce the length of the run to make it go a little faster. `threshold = ...` sets the maximum value of the multivariate Gelman Diagnostic used to assess convergence -- it defaults to 1.1, but we set it higher here for demonstrative purposes).

```{r}
invert.options$n.tries <- 1      # Number of attempts
invert.options$nchains <- 2      # Number of MCMC chains
invert.options$ngibbs <- 5000    # Number of iterations per chain
invert.options$burnin <- 1000    # Length of burnin period
invert.options$do.lsq.first <- TRUE # Initialize with results from a fast least-squares optimization algorithm
invert.options$threshold <- 1.3  # Maximum value for Gelman diagnostic
```

The full list of inversion options is below. See the documentation for `invert.auto` for a full description of each of these.

```{r}
names(invert.options)
```

Now, we load some test data (for _Acer rubrum_ leaves (`testspec_ACRU`, accessible through `data(testspec)`).

```{r}
data(testspec)
observed <- testspec_ACRU[,1]
plot(wl, observed, xlab="Wavelength", ylab="Reflectance", type='l')
```

To perform an inversion, just pass the observation matrix and the inversion settings into `invert.auto`. `quiet = TRUE` suppresses the progress bar, which is done here to clean up the knitted document. We also set `parallel = FALSE` here because this vignette is rebuilt on shared hardware where only one processor might be available; for runs on your own machine, you will probably want to set it TRUE to use all the cores on your machine.

```{r}
if(file.exists("inversion.output.rds")){
    inversion.output <- readRDS("inversion.output.rds")
} else {
    inversion.output <- invert.auto(observed = observed,
                                    invert.options = invert.options,
                                    quiet = TRUE,
                                    parallel = FALSE)
    saveRDS(inversion.output, "inversion.output.rds")
}
```

The output of `invert.auto` is a list of two objects: A list of summary statistics for each parameter and an `mcmc.list` object of the samples for diagnostic purposes or calculation of other summary statistics. 

```{r}
par(mfrow=c(2,2))
plot(inversion.output$samples, auto.layout=FALSE)

par(mfrow=c(1,1))
samples.mat <- as.matrix(inversion.output$samples)[-(2000:0),1:5]
colnames(samples.mat) <- params.prospect5
pairs(samples.mat, pch=".")

means <- unlist(inversion.output$results[grep("mu", names(inversion.output$results))])[1:5]
prospect.sim <- prospect(means, 5)[,1]  # reflectance

plot(wl, observed, type='l', col=1, xlab="wavelength (nm)", ylab="reflectance")
lines(wl, prospect.sim, type='l', col=2)
legend("topright", c("observed", "predicted"), lty=1, col=1:2)
```