---
title: "Tutorial: How to Use kwb.resilience"
author: "Andreas Matzinger"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7, 
  fig.height = 5
)
```

## Install the Package from GitHub

```{r eval = FALSE}
# install.packages("devtools")

remotes::install_github("kwb-r/kwb.resilience", dependencies = TRUE, build_vignettes = TRUE)
```

note: this will also install two handy KWB packages kwb.event and kwb.plot and their CRAN dependencies.

## Load the Package

```{r}
library("kwb.resilience")
```

## Define Helper Functions

We use the following plot function within this tutorial to look at our data:

```{r}
plot_dissolved_oxygen_S4 <- function(oxygen, ...) {
  plot(
    x = oxygen$timestamp, 
    y = oxygen$S4_red_Imp_Surface, 
    xlab = "Time", 
    ylab = "DO in mg/l",
    ...
  )
}
```

## Example Data

Let's have a look at some example data, which is included in kwb.resilience:

```{r}
# Show the first lines of the example data
head(oxygen)

# Plot the example data
plot_dissolved_oxygen_S4(oxygen)
```

## Use the Package

The following sections demonstrate the three main functions of the package:

* [resilience.severity](#severity)
* [resilience.events](#events)
* [resilience.summary](#summary)

### Severity and Resilience {#severity}

The function resilience.severity calculates severity as outlined in [Matzinger et al. (2018)](https://www.researchgate.net/publication/326040304_Quantitative_Beschreibung_der_Resilienz_urbaner_Wassersysteme).
The following example applies the function for scenario 4 of the test data set, using Pa = 2 mg/L as a lower threshold (for river fish) and Pmax = 0 mg/L as worst possible performance. Note: The function also works for upper thresholds.

```{r}
Sev.S4 <- resilience.severity(
  time_stamp = oxygen$timestamp, 
  Pt = oxygen$S4_red_Imp_Surface, 
  Pa = 2, 
  Pmax = 0
)

Sev.S4
```

The resilience index Res0 can then be calculated as 1-severity:

```{r}
Res0.S4 <- 1 - Sev.S4

Res0.S4
```


### Resilience Indices Separate for Each Event {#events}

A time series can contain several events when performance P(t) violates acceptable threshold Pa. The function resilience.events allows calculation of resilience indices separately for each event. Each line in the resulting data.frame shows one event. Again, let's make an example for scenario 4:

```{r}
# Define the acceptable performance
Pa <- 2

# Calculate Events
events.S4 <- resilience.events(
  time_stamp = oxygen$timestamp, 
  Pt = oxygen$S4_red_Imp_Surface, 
  Pa = Pa, 
  Pmax = 0, 
  evtSepTime = 6 * 60 * 60, 
  signalWidth = 15 * 60
)

events.S4
```

Let's plot the two failure events of scenario 4 (beginning and end is marked by vertical grey lines). The red line shows Pa = 2 mg/L. For the evaluation of resilience indices severity and Res0, only the area below Pa is integrated (in the case of an upper threshold it would be the area above Pa).

```{r}

# Plot Events
library(kwb.utils)

for (i in seq_len(nrow(events.S4))) {
  
  event <- events.S4[i, ]
  
  xlim <- kwb.utils::extendLimits(kwb.event::eventToXLim(event), 5)
  
  plot_dissolved_oxygen_S4(oxygen, xlim = xlim, ylim = c(0,6), 
                           main = paste("Event", i))
  
  abline(v = unlist(event[, c("tBeg", "tEnd")]), col = "grey")
  
  abline(h = Pa, col = "red")
  
}
```



### Resilience indices for entire time series {#summary}

The function resi.summary allows calculating resilience indices from [Matzinger et al. (2018)](https://www.researchgate.net/publication/326040304_Quantitative_Beschreibung_der_Resilienz_urbaner_Wassersysteme) for entire time series. The function combines complete integration of time series (e.g., for Res0) with event means (such as mean recovery time in per cent mean_trec_percent). In the case of resilience.summary all the time series (scenarios in test data) are evaluated. Each line of result shows result for one time series:


```{r}
resi.summary <- resilience.summary(
  time_stamp = oxygen$timestamp, 
  Pt = oxygen[, setdiff(names(oxygen), "timestamp")], 
  Pa = 2, 
  Pmax = 0, 
  evtSepTime = 6 * 60 * 60, 
  signalWidth = 15 * 60
)

resi.summary
```