euromomo: EuroMoMo algorithm

#' Generate outputs from the EuroMomo algorithm
#'
#' The function writes output graph files, tables into a textfile, etc.
#'
#' @param data \code{data.frame} generated by baseline function
#'
#' @section graphs:
#' graphs with crude mortality and z-scores
#' @section table:
#' table with cumulative mortality per season
#' @section data:
#' text file containing the original dataframe
#' @export
output <- function(data) {
  #
  # Graph: crude mortality
  #

  # Open connection to png file
  # Default size = 480 px, pointsize = 12
  filename <- paste0("Graph-Number_of_deaths-", groupOpts["label"], ".png")
  png(filename = file.path(week.dir, "output", filename), width = 800, height = 600, units = "px", pointsize = 12*600/480)

  # Create layout: one for the graph, one for the legend
  layout(mat = matrix(c(1, 2), nrow = 2, ncol = 1), heights = c(4, 1))
  # Graph
  par(mar = c(1.5, 2.5, 3, 2.5))
  plot.new()
  with(data, {
    plot.window(xlim = c(1, nrow(data)), ylim = range(cnb))
    idx <- which(WoDi == WoDi[nrow(data)])
    axis(side = 1, at = idx, labels = ISOweek[idx])
    axis(side = 2)
    box()
    title(main = paste("Number of deaths -", ISOweek[nrow(data)], "\n", getOption("euromomo")$Country, "-", groupOpts["label"]))
    # Add the graphs
    lines(x = 1:nrow(data), y = onb, col = "black")
    lines(x = 1:nrow(data), y = ifelse(onb != cnb, cnb, NA), col = "darkgreen")
    lines(x = 1:nrow(data), y = ifelse(cond == 1, onb, NA), col = "blue")
    lines(x = 1:nrow(data), y = pnb, col = "red")
    # Add prediction intervals
    # Always plot -2 and +2, and only plot prediction intervals that are smaller than max(onb)
    for (multiplier in c(-2, seq(from = 2, to = 20, by = 2))) {
      z <- (pnb^(2/3)+ multiplier*pv.pnb)^(3/2)
      if (is.element(multiplier, c(-2, 2)) | all(z < max(onb))) lines(x = 1:nrow(data), y = z, col = "orange")
    }
  })
  # Legend
  par(mar = rep(0, 4))
  plot.new()
  plot.window(xlim = c(0, 1), ylim = c(0, 1))
  legend(
    x = "center",
    xjust = 0.5, yjust = 0.5,
    legend = c("Known number of deaths", "Data used in model", "Corrected number deaths",
      "Baseline", "Prediction interval by 2 stdv"),
    lty = 1,
    col = c("black", "blue", "darkgreen", "red", "orange"),
    ncol = 2
  )

  # Close connection to png file
  dev.off()

  #
  # Graph: z-scores
  #

  # Open connection to png file
  # Default size = 480 px, pointsize = 12
  filename <- paste0("Graph-Zscore-", groupOpts["label"], ".png")
  png(filename = file.path(week.dir, "output", filename), width = 800, height = 600, units = "px", pointsize = 12*600/480)

  # Create layout: one for the graph, one for the legend
  layout(mat = matrix(c(1, 2), nrow = 2, ncol = 1), heights = c(4, 1))
  # Graph
  par(mar = c(1.5, 2.5, 3, 2.5))
  plot.new()
  with(data, {
    plot.window(xlim = c(1, nrow(data)), ylim = range(Zscore))
    idx <- which(WoDi == WoDi[nrow(data)])
    axis(side = 1, at = idx, labels = ISOweek[idx])
    axis(side = 2, at = seq(from = -20, to = 20, by = 2))
    abline(h =  seq(from = -20, to = 20, by = 2), col = "orange")
    abline(h = 0, col = "red")
    box()
    title(main = paste("Z-score -", ISOweek[nrow(data)], "\n", getOption("euromomo")$Country, "-", groupOpts["label"]))
    # Add the graphs
    lines(x = 1:nrow(data), y = Zscore, col = "black")
    lines(x = 1:nrow(data), y = ifelse(cond == 1, Zscore, NA), col = "blue")
  })
  # Legend
  par(mar = rep(0, 4))
  plot.new()
  plot.window(xlim = c(0, 1), ylim = c(0, 1))
  legend(
    x = "center",
    xjust = 0.5, yjust = 0.5,
    legend = c("Z-score", "Data used in model", "Baseline", "Standard deviations by 2"),
    lty = 1,
    col = c("black", "blue", "red", "orange"),
    ncol = 2
  )

  # Close connection to png file
  dev.off()

  #
  # Table with cumulative mortality per season
  #

  # Define variable that defines winter and summer season
  data <- within(data, season <- ifelse(WoDi >= 40 | WoDi <= 20, paste0("Winter-", YoDi-(WoDi <= 20)), paste0("Summer-", YoDi)))

  # Split data according to season -> generates a list of dataframes
  data.split <- with(data, split(data, f = season))

  # Apply cumulative moratality function to each dataframe in the list data.split
  tmp <- lapply(data.split, function(x) {
    with(x, data.frame(
      Start = ISOweek[1],
      End = ISOweek[nrow(x)],
      Nweeks = nrow(x),
      CumDeaths = round(sum(cnb), digits = 2),
      CumExcess = round(sum(excess), digits = 2)))
  })

  # And apppend them togeter into one dataframe
  cummort.data <- do.call(what = "rbind", args = tmp)

  # Save cummort.data in a text file
  filename <- paste0("Table-Cumulative_mortality-", groupOpts["label"], ".txt")
  capture.output(print(cummort.data), file = file.path(week.dir, "output", filename))

  #
  # Text file containing the original dataframe
  #

  filename <- paste0("Data-", groupOpts["label"], ".txt")
  write.table(x = data, file = file.path(week.dir, "output", filename),
    quote = FALSE, sep = ";", row.names = FALSE)
}

#' Write final output to HUB compatible files.
#'
#' @param final Final data.frame
#' @param dLastFullWeek Last final week.
#' @return Nothing
#' @export
#'
writeHUBFiles <- function(final, dLastFullWeek) {
  #Restrict to EurMOMO hub pre-specified groups.
  final4hub.restricted <- final4hub.complete <- subset(final, group.name %in% paste("momodefault",1:5,sep=""))
  #Restrict to "legal" columns.
  final4hub.restricted[,c("nb","onb","cnb","pnb")] <- NA

  #Build up the file name.
  countryStr <- getOption("euromomo")[["Country"]]
  ISOweek <- ISOweek(momoFile$dLastFullWeek)
  write.table(final4hub.complete, file = file.path(week.dir, "data", filename=paste(countryStr,"-complete-",ISOweek,".txt",sep="")), quote = FALSE, sep = ";", row.names = FALSE)
  write.table(final4hub.restricted, file = file.path(week.dir, "data", filename=paste(countryStr,"-restricted-",ISOweek,".txt",sep="")), quote = FALSE, sep = ";", row.names = FALSE)

  #Small check.
  table(final4hub.complete$group.name)
  head(final4hub.restricted)

  invisible()
}