---
title: "02 Getting started with `SpaDES`"
author: "Eliot J. B. McIntire"
date: "`r format(Sys.Date(), format='%B %d %Y')`"
output:
  pdf_document:
  number_sections: yes
  toc: yes
html_document:
  self_contained: no
  toc: yes
vignette: >
  %\VignetteEngine{knitr::knitr}
  %\VignetteIndexEntry{02 Getting started with SpaDES}
  %\VignetteDepends{SpaDES}
  %\VignetteKeyword{tutorial}
  %\usepackage[utf8]{inputenc}
---

# Getting started with `SpaDES`

- module organization/overview; module development philosophy

-----

# Sample model

```{r using-SpaDES, eval=TRUE, echo=TRUE, fig=TRUE, cache=TRUE}
library(SpaDES)

times <- list(start = 0.0, end = 20.0)
parameters <- list(
  .globals = list(stackName = "landscape", burnStats = "nPixelsBurned"),
  .progress = list(NA),
  randomLandscapes = list(nx = 100L, ny = 100L, inRAM = TRUE),
  fireSpread = list(nFires = 10L, spreadprob = 0.225, its = 1e6, persistprob = 0,
                    returnInterval = 10, startTime = 0,
                    .plotInitialTime = 0.0, .plotInterval = 10),
  caribouMovement = list(N = 100L, moveInterval = 1,
                         .plotInitialTime = 1, .plotInterval = 1)
)
modules <- list("randomLandscapes", "fireSpread", "caribouMovement")
path <- list(modulePath = system.file("sampleModules", package = "SpaDES.core"),
             outputPath = file.path(tempdir(), "simOutputs"))

mySim <- simInit(times = times, params = parameters, modules = modules, path = path)

if (interactive()) dev()
spades(mySim)
```

## fire spread
```{r fire, eval=TRUE, echo=TRUE, fig=TRUE, cache=TRUE, message=FALSE}
library(ggplot2)
library(RColorBrewer)
nFires <- 10L
landscape[["Fires"]] <- spread(landscape[["percentPine"]],
                               loci = sample(1L:ncell(landscape), nFires),
                               spreadProb = landscape[["percentPine"]] / 
                                 (maxValue(landscape[["percentPine"]]) * 5) + 0.1,
                               persistance = 0, mask = NULL, maxSize = 1e8,
                               directions = 8, iterations = 1e6, plot.it = FALSE,
                               mapID = TRUE)
setColors(landscape$Fires) <- brewer.pal(8, "Reds")[5:8]

Plot(landscape[["Fires"]], new = TRUE)
fireSizes <- getValues(landscape[["Fires"]])
ggHist <- qplot(fireSizes[fireSizes > 0])
Plot(ggHist)
```

```{r fire-overlaid, eval=TRUE, echo=TRUE, fig=TRUE}
# Show the burning more strongly over abundant pine
percentPine <- landscape$percentPine
Plot(percentPine, new = TRUE)
Plot(landscape[["Fires"]], addTo = "percentPine", legend = FALSE, title = FALSE)
```

We can see that the fires tend to be in the Pines because we made it that way, using an arbitrary weighting with pine abundance:

```{r fire-impacts, eval=TRUE, echo=TRUE, fig=FALSE}
# Show the burning more strongly over abundant pine
fire <- reclassify(landscape[["Fires"]], rcl = cbind(0:1, c(0,ncell(landscape)), 0:1))
pine <- reclassify(landscape[["percentPine"]], rcl = cbind(0:9*10, 1:10*10, 0:9))
PineByFire <- crosstab(fire, pine, long=TRUE)
colnames(PineByFire) <- c("fire", "pine", "freq")
PineByFire$pine <- as.numeric(as.character(PineByFire$pine))
summary(glm(freq ~ fire*pine, data=PineByFire, family="poisson"))
```

Sure enough, there are more fires as the abundance of pine goes up, as seen by the positive interaction term (the negative `fire1` term means that there are more pixels without fires than with fires).

**Impact some of the forest**

```{r fire-impacts-maps, eval=TRUE, echo=TRUE, fig=TRUE}
landscape[["forestAge"]][landscape[["Fires"]] > 0] <- 0
landscape[["habitatQuality"]][landscape[["Fires"]] > 0] <- 0.1
landscape[["percentPine"]][landscape[["Fires"]] > 0] <- 0
Plot(landscape, new = TRUE)
```

To model point agents (mobile *e.g.*, animals, or non-mobile, *e.g.*, plants), use a `SpatialPointsDataFrame` containing additional columns for storing agents' previous `n` positions.

```{r mobile-point-agent, echo=TRUE, eval=TRUE}
N <- 10L # number of agents

# caribou data vectors
IDs <- letters[1L:N]
sex <- sample(c("female", "male"), N, replace = TRUE)
age <- round(stats::rnorm(N, mean = 8, sd = 3))
x1 <- stats::runif(N, -50, 50) # previous X location
y1 <- stats::runif(N, -50, 50) # previous Y location

# caribou (current) coordinates
x0 <- stats::rnorm(N, x1, 5)
y0 <- stats::rnorm(N, y1, 5)

# create the caribou agent object
caribou <- data.frame(x1, y1, sex, age)
coordinates(caribou) <- cbind(x = x0, y = y0)
row.names(caribou) <- IDs
```

Using a simple landscape-dependent correlated random walk, we simulate the movement of caribou across a heterogeneous landscape. Because we had just had fires, and we assume that fires have a detrimental effect on animal movement, we can see the long steps taken in the new, low quality, post-burn sections of the landscape.

```{r agent-crw-trajectory, eval=TRUE, echo=TRUE, fig=TRUE}
library(sp)
if (interactive()) dev()
Plot(landscape[["habitatQuality"]], new = TRUE)

for (t in 1:30) {
  #crop any caribou that went off maps
  caribou <<- crop(caribou,landscape)
  caribouArrows <- makeLines(from = SpatialPoints(cbind(x = caribou$x1, y = caribou$y1)),
            to = caribou)
  Plot(caribouArrows, length = 0.04, addTo = "landscape$habitatQuality")

  # find out what pixels the individuals are on now
  ex <- landscape[["habitatQuality"]][caribou]

  #step length is a function of current cell's landscape quality
  sl <- 0.25/ex
  ln <- rlnorm(length(ex), sl, 0.02) # log normal step length
  sd <- 30

  caribou <<- crw(caribou, stepLength = ln, stddev = sd, lonlat = FALSE)
}
```

Below is a greatly simplified example, using the sample `randomLandscapes` and `fireSpread` modules. *NB only two parameters are varied; no outputs are saved; and the analyses done here are kept simple for illustrative purposes. This will take a while to run!*

```{r multiple-simulations, echo=TRUE, eval=FALSE}
###
### WARNING: this can take a while to run, especially for large mapSizes.
###
rasterOptions(maxmemory = 1e9)

# list all parameter values with which to run simulations
parameters <- list(mapSize = round(sqrt(c(1e4, 1e5, 1e6, 1e7, 1e8))),
                   pSpread = seq(0.05, 0.25, 0.05))

# create data.frame with all parameter combinations
paramsdf <- expand.grid(parameters)

# outputs
nPixelsBurned <- numeric()
meanPixelsBurned <- cbind(paramsdf, pmean = NA, psd = NA)

set.seed(42)
for (i in 1:nrow(paramsdf)) {
  # initialize each simulation with a param combo from paramsdf
  mySim <- with(paramsdf,
                simInit(times = list(start = 0.0, end = 20.0),
                        params = list(
                          .progress = list(type = NA, interval = NA),
                          .globals = list(stackName = "landscape", burnStats = "nPixelsBurned"),
                          randomLandscapes = list(nx = mapSize[i], ny = mapSize[i],
                                                inRAM = TRUE),
                          fireSpread = list(nFires = 1000, spreadprob = pSpread[i],
                                          persistprob = 0, its = 1e6,
                                          returnInterval = 10, startTime = 0)
                          ),
                        modules = list("randomLandscapes", "fireSpread"),
                        path = system.file("sampleModules", package = "SpaDES.core")))
  mySim <- spades(mySim)

  # collect stats for each run
  proportionBurned <- with(paramsdf, nPixelsBurned / (mapSize[i]^2))
  meanPixelsBurned[i, "pmean"] <- mean(proportionBurned)
  meanPixelsBurned[i, "psd"] <- sd(proportionBurned)

  # cleanup between runs
  rm(landscape, mySim, nPixelsBurned)
  for (j in 1:10) gc()
}

# overall statistics
pch <- c(21:25)
col <- brewer.pal(5, "Set1")

with(meanPixelsBurned, plot(pmean ~ pSpread, xlab = "Spread probability",
                            ylab = "Mean proportion of pixels burned",
                            ylim = c(0,1), pch = pch, cex = 1.5, col = col))
with(parameters, legend("topleft", legend = formatC(mapSize^2, digits = 0),
                        pch = pch, col = col, cex = 1.2))
```

```{r multiple-simulations-outputs, echo=FALSE, eval=FALSE}
# this is included as the output from the previous chunk,
# so you don't need to wait for the prev chunk to run
meanPixelsBurned <- structure(list(mapSize = c(100, 316, 1000, 3162, 10000, 100,
316, 1000, 3162, 10000, 100, 316, 1000, 3162, 10000, 100, 316,
1000, 3162, 10000, 100, 316, 1000, 3162, 10000), pSpread = c(0.05,
0.05, 0.05, 0.05, 0.05, 0.1, 0.1, 0.1, 0.1, 0.1, 0.15, 0.15,
0.15, 0.15, 0.15, 0.2, 0.2, 0.2, 0.2, 0.2, 0.25, 0.25, 0.25,
0.25, 0.25), pmean = c(0.1494, 0.0160831597500401, 0.00157566666666667,
0.000160061440122219, 1.57933333333333e-05, 0.241933333333333,
0.0301600972066442, 0.00304333333333333, 0.000307687363234317,
3.14933333333333e-05, 0.4268, 0.0709388185654009, 0.00855766666666667,
0.000826111731886786, 8.185e-05, 0.673466666666667, 0.271898200074774,
0.0522476666666667, 0.0057268723054435, 0.000569923333333333,
0.8551, 0.704143967312931, 0.758393333333333, 0.748152608931462,
0.755049233333333), psd = c(0.000529150262212918, 0.00034444366932222,
2.11266025033211e-05, 5.52032923607854e-06, 2.80416357107309e-07,
0.00248461935381122, 0.00117238603350013, 8.3500499000505e-05,
5.02117387430545e-06, 6.50640709864774e-08, 0.0103764155660806,
0.00292733584237083, 0.000185057648675577, 2.23794372722204e-05,
3.33725935462019e-06, 0.00638931399551884, 0.0181866853325732,
0.00354059674819561, 6.08251802692677e-05, 3.54532908674686e-05,
0.00927739187487513, 0.00349463109959475, 0.000899981296101935,
0.00115857888357641, 0.000319090208300621)), .Names = c("mapSize",
"pSpread", "pmean", "psd"), row.names = c(NA, -25L), class = "data.frame")
```

![](../inst/figs/plot-area-burned.pdf "Mean proportion of pixels burned for maps of various sizes and fire spread probabilities.")

---

### loading files (overlaps with modules vignette?)

Bringing data into `R` is easy, and can be done using any of the built in data import tools. To facilitate this, we have provided additional functionality to easily load maps or data from files via the load module. To automatically import a list of files, simply provide it as a parameter named \texttt{.load$fileList} when initializing the simulation. See \texttt{?loadFiles} and the modules vignette for more information on the load module.

```{r load-landscape-maps, echo=TRUE, eval=TRUE, cache=TRUE, fig=TRUE}
### Example: loading habitat maps

# use all built-in maps from the SpaDES package
pathToMaps <- file.path(find.package("quickPlot", quiet = FALSE), "maps")
fileList <- data.frame(files = dir(pathToMaps, full.names = TRUE, pattern = "tif"),
                       functions = "rasterToMemory", packages = "SpaDES.core",
                       stringsAsFactors = FALSE)

# this list can be passed to simInit() as an entry in the parameter list
mySim <- simInit(times = list(start = 0.0, end = 10),
                 params = list(
                   .load = list(fileList = fileList),
                   .progress = list(NA),
                   .globals = list(stackName = "landscape", burnStats = "nPixelsBurned")
                 ),
                 modules = list("fireSpread"),
                 path = system.file("sampleModules", package = "SpaDES.core"))

# put layers into a single stack for convenience
landscape <- stack(DEM, forestCover, forestAge, habitatQuality, percentPine)

spades(mySim)
```

# Debugging simulations

https://github.com/PredictiveEcology/SpaDES/wiki/Debugging
