One box time series and twa values

- one_box() creates decline time series from mkinfit objects or simply
  from a half-life
- sawtooth() generates sawtooth curves for arbitrary application
  patterns and decline models
- twa() calculates moving window averages
- max_twa() gives their maxima and
- plot.one_box() can plot series generated by one_box() or sawtooth(),
  optionally adding a greay rectangle to illustrate the maximum moving
  window time weighted average

9 files changed, 423 insertions, 3 deletions

diff --git a/ChangeLog b/ChangeLog
index e78b175..fd5074e 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,20 @@
+commit bba2cf3a70849ba86f37520d3e909cf1c706f416
+Author: Johannes Ranke <jranke@uni-bremen.de>
+Date:   2016-12-22 11:06:57 +0100
+
+    Fix reading in times from .out files
+    
+    The code from the previous commit was broken. Also, the time
+    zone for the times that are read is now wet to 'UTC', in order to
+    avoid setting different time zones due to daylight savings, which
+    introduces artificial one-hour offsets on changeover days.
+
+commit 0af7c7b8c34067fc4756929925239c329b28ed32
+Author: Johannes Ranke <jranke@uni-bremen.de>
+Date:   2016-12-14 18:26:14 +0100
+
+    Changelog update and roxygen run
+
 commit 5a04ad3061c1484b45703e44149f49ec97cfbf15
 Author: Johannes Ranke <jranke@uni-bremen.de>
 Date:   2016-12-14 16:52:14 +0100
diff --git a/DESCRIPTION b/DESCRIPTION
index 7207ca4..f92e6da 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: pfm
 Type: Package
 Title: Utilities for Pesticide Fate Modelling
-Version: 0.3-9
-Date: 2016-12-22
+Version: 0.4-1
+Date: 2017-01-18
 Authors@R: person("Johannes Ranke", email = "jranke@uni-bremen.de", 
                   role = c("aut", "cre", "cph"))
 Description: Utilities for simple calculations of predicted environmental
@@ -17,7 +17,8 @@ Imports:
     methods
 Suggests:
     testthat,
-    chents
+    chents,
+    magrittr
 License: GPL
 LazyLoad: yes
 LazyData: yes
diff --git a/NAMESPACE b/NAMESPACE
index 3777574..af0eefa 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,8 +2,13 @@
 
 S3method(GUS,chent)
 S3method(GUS,numeric)
+S3method(max_twa,one_box)
+S3method(one_box,mkinfit)
+S3method(one_box,numeric)
 S3method(plot,TOXSWA_cwa)
+S3method(plot,one_box)
 S3method(print,GUS_result)
+S3method(twa,one_box)
 export(GUS)
 export(PEC_soil)
 export(PEC_sw_drainage_UK)
@@ -14,15 +19,25 @@ export(SSLRC_mobility_classification)
 export(TOXSWA_cwa)
 export(endpoint)
 export(geomean)
+export(max_twa)
+export(one_box)
 export(pfm_degradation)
 export(read.TOXSWA_cwa)
+export(sawtooth)
 export(soil_DT50)
 export(soil_Kfoc)
 export(soil_N)
 export(soil_sorption)
+export(twa)
 import(graphics)
 import(mkin)
 importFrom(R6,R6Class)
 importFrom(methods,is)
+importFrom(mkin,mkinpredict)
 importFrom(readr,fwf_empty)
 importFrom(readr,read_fwf)
+importFrom(stats,filter)
+importFrom(stats,frequency)
+importFrom(stats,plot.ts)
+importFrom(stats,time)
+importFrom(stats,ts)
diff --git a/R/twa.R b/R/twa.R
new file mode 100644
index 0000000..3865fd7
--- /dev/null
+++ b/R/twa.R
@@ -0,0 +1,220 @@
+# Copyright (C) 2017  Johannes Ranke
+
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Create a time series of decline data
+#'
+#' The time series starts with the amount specified for the first application.
+#' This does not create objects of type \code{\link{ts}}.
+#'
+#' @param x When numeric, this is the half-life to be used for an exponential
+#' decline. If x is an mkinfit object, the decline is calculated from this object
+#' @param t_end End of the time series
+#' @param res Resolution of the time series
+#' @param ... Further arguments passed to methods
+#' @importFrom stats filter frequency time ts
+#' @export
+#' @examples
+#' # Only use a half-life
+#' pred_0 <- one_box(10)
+#' plot(pred_0)
+#'
+#' # Use a fitted mkinfit model
+#' require(mkin)
+#' fit <- mkinfit("FOMC", FOCUS_2006_C)
+#' pred_1 <- one_box(fit)
+#' plot(pred_1)
+#'
+#' # Use a model with more than one observed variable
+#' m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
+#' fit_2 <- mkinfit(m_2, FOCUS_2006_D)
+#' pred_2 <- one_box(fit_2)
+#' plot(pred_2)
+one_box <- function(x,
+  t_end = 100, res = 0.01, ...)
+{
+  UseMethod("one_box")
+}
+
+#' @rdname one_box
+#' @export
+one_box.numeric <- function(x,
+  t_end = 100, res = 0.01, ...)
+{
+  half_life = x
+  k = log(2)/half_life
+  t_out <- seq(0, t_end, by = res)
+  raw <- matrix(exp( - k * t_out), ncol = 1)
+  dimnames(raw) <- list(NULL, "parent")
+  result <- ts(raw, 0, t_end, frequency = 1/res)
+  class(result) <- c("one_box", "ts")
+  return(result)
+}
+
+#' @rdname one_box
+#' @importFrom mkin mkinpredict
+#' @export
+one_box.mkinfit <- function(x, t_end = 100, res = 0.01, ...) {
+  fit <- x
+  t_out = seq(0, t_end, by = res)
+  odeini <- c(1, rep(0, length(fit$mkinmod$spec) - 1))
+  names(odeini) <- names(fit$mkinmod$spec)
+  if (length(fit$mkinmod$spec) == 1) solution_type = "analytical"
+  else solution_type = "deSolve"
+
+  tmp <- mkinpredict(fit$mkinmod, odeparms = fit$bparms.ode, odeini = odeini,
+    outtimes = t_out, solution_type = solution_type)[-1]
+  result <- ts(tmp, 0, t_end, frequency = 1/res)
+  class(result) <- c("one_box", "ts")
+  return(result)
+}
+
+#' Plot time series of decline data
+#' 
+#' @param x The object of type \code{\link{one_box}} to be plotted
+#' @param xlim Limits for the x axis
+#' @param ylim Limits for the y axis
+#' @param xlab Label for the x axis
+#' @param ylab Label for the y axis
+#' @param max_twa If a numeric value is given, the maximum time weighted
+#'   average concentration(s) is/are shown in the graph.
+#' @param max_twa_var Variable for which the maximum time weighted average should
+#'   be shown if max_twa is not NULL.
+#' @param ... Further arguments passed to methods
+#' @importFrom stats plot.ts
+#' @export
+#' @examples
+#' plot(sawtooth(one_box(10), 3, 7), max_twa = 21)
+plot.one_box <- function(x, 
+                         xlim = range(time(x)), ylim = c(0, max(x)),
+                         xlab = "Time", ylab = "Fraction of initial",
+                         max_twa = NULL, max_twa_var = dimnames(x)[[2]][1], ...)
+{
+  obs_vars <- dimnames(x)[[2]]
+  plot.ts(x, plot.type = "single", xlab = xlab, ylab = ylab,
+          lty = 1:length(obs_vars), col = 1:length(obs_vars),
+          las = 1, xlim = xlim, ylim = ylim)
+  if (!is.null(max_twa)) {
+    x_twa <- max_twa(x, window = max_twa)
+    value <- x_twa$max[max_twa_var]
+    rect(x_twa$window_start[max_twa_var], 0, 
+         x_twa$window_end[max_twa_var], value, col = "grey")
+    text(x_twa$window_end[max_twa_var], value, paste("Maximum:", signif(value, 3)), pos = 4)
+    # Plot a second time to cover the grey rectangle
+    matlines(time(x), as.matrix(x), lty = 1:length(obs_vars), col = 1:length(obs_vars))
+  }
+}
+
+#' Create decline time series for multiple applications
+#'
+#' If the number of application cycles \code{n} is greater than 1, the
+#' application pattern specified in \code{applications} is repeated \code{n}
+#' times, with an interval \code{i}.
+#' @param x A \code{\link{one_box}} object
+#' @param n The number of applications. If \code{applications} is specified, \code{n} is not used
+#' @param i The interval between applications. If \code{applications} is specified, \code{i} 
+#'   is not used
+#' @param applications A data frame holding the application times in the first column and
+#'   the corresponding amounts applied in the second column for each application cycle.
+#'   If \code{n} is one, the application pattern specified here is used only once.
+#' @export
+#' @examples
+#' applications = data.frame(time = seq(0, 14, by = 7), amount = c(1, 2, 3))
+#' pred <- one_box(10)
+#' plot(sawtooth(pred, applications = applications))
+#'
+#' m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
+#' fit_2 <- mkinfit(m_2, FOCUS_2006_D)
+#' pred_2 <- one_box(fit_2)
+#' pred_2_saw <- sawtooth(pred_2, 2, 7)
+#' plot(pred_2_saw, max_twa = 21, max_twa_var = "m1")
+#'
+#' max_twa(pred_2_saw)
+sawtooth <- function(x, n = 1, i = 365, 
+                     applications = data.frame(time = seq(0, 0 + n * i, length.out = n), 
+                                               amount = 1))
+{
+  n_obs = ncol(as.matrix(x))
+  t_end = max(time(x))
+  freq = frequency(x)
+  empty <- ts(matrix(0, nrow = t_end * freq, ncol = n_obs), 0, t_end, freq)
+  result <- empty
+  for (i_app in 1:nrow(applications)) {
+    t_app <- applications[i_app, "time"]
+    amount_app <- applications[i_app, "amount"]
+    if (t_app == 0) {
+      result <- result + x * amount_app
+    } else {
+      lag_phase <- as.matrix(empty)[1:(t_app * freq), , drop = FALSE]
+      app_phase <- amount_app * as.matrix(x)[1:((t_end - t_app) * freq + 1), , drop = FALSE]
+      app_ts <- ts(rbind(lag_phase, app_phase), 0, t_end, frequency = freq)
+      result <- result + app_ts
+    }
+  }
+  class(result) = c("one_box", "ts")
+  dimnames(result) <- dimnames(x)
+  return(result)
+}
+
+#' Calculate a time weighted average concentration
+#'
+#' The moving average is built only using the values in the past, so
+#' the earliest possible time for the maximum in the time series returned
+#' is after one window has passed.
+#'
+#' @param x An object of type \code{\link{one_box}}
+#' @param window The size of the moving window
+#' @seealso max_twa
+#' @export
+#' @examples
+#' pred <- sawtooth(one_box(10),
+#'   applications = data.frame(time = c(0, 7), amount = c(1, 1)))
+#' max_twa(pred)
+twa <- function(x, window = 21) UseMethod("twa")
+
+#' @rdname twa
+#' @export
+twa.one_box <- function(x, window = 21)
+{
+  resolution = 1/frequency(x)
+  n_filter = window/resolution
+  result = filter(x, rep(1/n_filter, n_filter), method = "convolution", sides = 1)
+  class(result) = c("one_box", "ts")
+  dimnames(result) <- dimnames(x)
+  return(result)
+}
+
+#' The maximum time weighted average concentration for a moving window
+#' 
+#' @seealso twa
+#' @inheritParams twa
+#' @export
+max_twa <- function(x, window = 21) UseMethod("max_twa")
+
+#' @export
+max_twa.one_box <- function(x, window = 21)
+{
+  freq = frequency(x)
+
+  twa_ts <- twa(x, window = window)
+  window_end <- apply(twa_ts, 2, which.max) / freq
+  result <- list()
+  result$max <- apply(twa_ts, 2, max, na.rm = TRUE)
+  result$window_start <- window_end - window
+  result$window_end <- window_end
+  return(result)
+}
diff --git a/man/max_twa.Rd b/man/max_twa.Rd
new file mode 100644
index 0000000..13e7bbf
--- /dev/null
+++ b/man/max_twa.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/twa.R
+\name{max_twa}
+\alias{max_twa}
+\title{The maximum time weighted average concentration for a moving window}
+\usage{
+max_twa(x, window = 21)
+}
+\arguments{
+\item{x}{An object of type \code{\link{one_box}}}
+
+\item{window}{The size of the moving window}
+}
+\description{
+The maximum time weighted average concentration for a moving window
+}
+\seealso{
+twa
+}
diff --git a/man/one_box.Rd b/man/one_box.Rd
new file mode 100644
index 0000000..86440ce
--- /dev/null
+++ b/man/one_box.Rd
@@ -0,0 +1,45 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/twa.R
+\name{one_box}
+\alias{one_box}
+\alias{one_box.numeric}
+\alias{one_box.mkinfit}
+\title{Create a time series of decline data}
+\usage{
+one_box(x, t_end = 100, res = 0.01, ...)
+
+\method{one_box}{numeric}(x, t_end = 100, res = 0.01, ...)
+
+\method{one_box}{mkinfit}(x, t_end = 100, res = 0.01, ...)
+}
+\arguments{
+\item{x}{When numeric, this is the half-life to be used for an exponential
+decline. If x is an mkinfit object, the decline is calculated from this object}
+
+\item{t_end}{End of the time series}
+
+\item{res}{Resolution of the time series}
+
+\item{...}{Further arguments passed to methods}
+}
+\description{
+The time series starts with the amount specified for the first application.
+This does not create objects of type \code{\link{ts}}.
+}
+\examples{
+# Only use a half-life
+pred_0 <- one_box(10)
+plot(pred_0)
+
+# Use a fitted mkinfit model
+require(mkin)
+fit <- mkinfit("FOMC", FOCUS_2006_C)
+pred_1 <- one_box(fit)
+plot(pred_1)
+
+# Use a model with more than one observed variable
+m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
+fit_2 <- mkinfit(m_2, FOCUS_2006_D)
+pred_2 <- one_box(fit_2)
+plot(pred_2)
+}
diff --git a/man/plot.one_box.Rd b/man/plot.one_box.Rd
new file mode 100644
index 0000000..7d30bd2
--- /dev/null
+++ b/man/plot.one_box.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/twa.R
+\name{plot.one_box}
+\alias{plot.one_box}
+\title{Plot time series of decline data}
+\usage{
+\method{plot}{one_box}(x, xlim = range(time(x)), ylim = c(0, max(x)),
+  xlab = "Time", ylab = "Fraction of initial", max_twa = NULL,
+  max_twa_var = dimnames(x)[[2]][1], ...)
+}
+\arguments{
+\item{x}{The object of type \code{\link{one_box}} to be plotted}
+
+\item{xlim}{Limits for the x axis}
+
+\item{ylim}{Limits for the y axis}
+
+\item{xlab}{Label for the x axis}
+
+\item{ylab}{Label for the y axis}
+
+\item{max_twa}{If a numeric value is given, the maximum time weighted
+average concentration(s) is/are shown in the graph.}
+
+\item{max_twa_var}{Variable for which the maximum time weighted average should
+be shown if max_twa is not NULL.}
+
+\item{...}{Further arguments passed to methods}
+}
+\description{
+Plot time series of decline data
+}
+\examples{
+plot(sawtooth(one_box(10), 3, 7), max_twa = 21)
+}
diff --git a/man/sawtooth.Rd b/man/sawtooth.Rd
new file mode 100644
index 0000000..4576c8b
--- /dev/null
+++ b/man/sawtooth.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/twa.R
+\name{sawtooth}
+\alias{sawtooth}
+\title{Create decline time series for multiple applications}
+\usage{
+sawtooth(x, n = 1, i = 365, applications = data.frame(time = seq(0, 0 + n
+  * i, length.out = n), amount = 1))
+}
+\arguments{
+\item{x}{A \code{\link{one_box}} object}
+
+\item{n}{The number of applications. If \code{applications} is specified, \code{n} is not used}
+
+\item{i}{The interval between applications. If \code{applications} is specified, \code{i} 
+is not used}
+
+\item{applications}{A data frame holding the application times in the first column and
+the corresponding amounts applied in the second column for each application cycle.
+If \code{n} is one, the application pattern specified here is used only once.}
+}
+\description{
+If the number of application cycles \code{n} is greater than 1, the
+application pattern specified in \code{applications} is repeated \code{n}
+times, with an interval \code{i}.
+}
+\examples{
+applications = data.frame(time = seq(0, 14, by = 7), amount = c(1, 2, 3))
+pred <- one_box(10)
+plot(sawtooth(pred, applications = applications))
+
+m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
+fit_2 <- mkinfit(m_2, FOCUS_2006_D)
+pred_2 <- one_box(fit_2)
+pred_2_saw <- sawtooth(pred_2, 2, 7)
+plot(pred_2_saw, max_twa = 21, max_twa_var = "m1")
+
+max_twa(pred_2_saw)
+}
diff --git a/man/twa.Rd b/man/twa.Rd
new file mode 100644
index 0000000..a6afa87
--- /dev/null
+++ b/man/twa.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/twa.R
+\name{twa}
+\alias{twa}
+\alias{twa.one_box}
+\title{Calculate a time weighted average concentration}
+\usage{
+twa(x, window = 21)
+
+\method{twa}{one_box}(x, window = 21)
+}
+\arguments{
+\item{x}{An object of type \code{\link{one_box}}}
+
+\item{window}{The size of the moving window}
+}
+\description{
+The moving average is built only using the values in the past, so
+the earliest possible time for the maximum in the time series returned
+is after one window has passed.
+}
+\examples{
+pred <- sawtooth(one_box(10),
+  applications = data.frame(time = c(0, 7), amount = c(1, 1)))
+max_twa(pred)
+}
+\seealso{
+max_twa
+}