# Copyright (C) 2019 Johannes Ranke
# Contact: jranke@uni-bremen.de

# This file is part of the R package mkin

# mkin 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/>

nafta <- function(ds, title = NA, quiet = FALSE, ...) {
  if (length(levels(ds$name)) > 1) {
    stop("The NAFTA procedure is only defined for decline data for a single compound")
  }
  n <- nrow(subset(ds, !is.na(value)))
  models <- c("SFO", "IORE", "DFOP")

  result <- list(title = title, data = ds)
  result$mmkin <- mmkin(models, list(ds), quiet = TRUE, ...)

  distimes <- lapply(result$mmkin, function(x) as.numeric(endpoints(x)$distimes["parent", ]))

  result$distimes <- matrix(NA, nrow = 3, ncol = 3,
    dimnames = list(models, c("DT50", "DT90", "DT50_rep")))
  result$distimes["SFO", ] <- distimes[[1]][c(1, 2, 1)]
  result$distimes["IORE", ] <- distimes[[2]][c(1, 2, 3)]
  result$distimes["DFOP", ] <- distimes[[3]][c(1, 2, 4)]

  # Get parameters with statistics
  result$parameters <- lapply(result$mmkin, function(x) {
    summary(x)$bpar[, c(1, 4:6)]
  })
  names(result$parameters) <- models

  # Compare the sum of squared residuals (SSR) to the upper bound of the
  # confidence region of the SSR for the IORE model
  result$S <- sapply(result$mmkin, function(x) sum(x$data$residual^2))
  names(result$S) <- c("SFO", "IORE", "DFOP")
  # Equation (3) on p. 3
  p <- 3
  result$S["IORE"]
  result$S_c <- result$S[["IORE"]] * (1 + p/(n - p) * qf(0.5, p, n - p))

  result$t_rep <- .evaluate_nafta_results(result$S, result$S_c,
    result$distimes, quiet = quiet)

  class(result) <- "nafta"
  return(result)
}

plot.nafta <- function(x, legend = FALSE, main = "auto", ...) {
  if (main == "auto") {
    if (is.na(x$title)) main = ""
    else main = x$title
  }
  plot(x$mmkin, ..., legend = legend, main = main)
}

print.nafta <- function(x, quiet = TRUE, digits = 3, ...) {
  cat("Sums of squares:\n")
  print(x$S)
  cat("\nCritical sum of squares for checking the SFO model:\n")
  print(x$S_c)
  cat("\nParameters:\n")
  print(x$parameters, digits = digits)
  t_rep <- .evaluate_nafta_results(x$S, x$S_c, x$distimes, quiet = quiet)
  cat("\nDTx values:\n")
  print(signif(x$distimes, digits = digits))
  cat("\nRepresentative half-life:\n")
  print(round(t_rep, 2))
}

.evaluate_nafta_results <- function(S, S_c, distimes, quiet = FALSE) {
  t_SFO <- distimes["IORE", "DT50"]
  t_IORE <- distimes["IORE", "DT50_rep"]
  t_DFOP2 <- distimes["DFOP", "DT50_rep"]

  if (S["SFO"] < S_c) {
    if (!quiet) {
      message("S_SFO is lower than the critical value S_c, use the SFO model")
    }
    t_rep <- t_SFO
  } else {
    if (!quiet) {
      message("The SFO model is rejected as S_SFO is equal or higher than the critical value S_c")
    }
    if (t_IORE < t_DFOP2) {
      if (!quiet) {
        message("The half-life obtained from the IORE model may be used")
      }
      t_rep <- t_IORE
    } else {
      if (!quiet) {
        message("The representative half-life of the IORE model is longer than the one corresponding")
        message("to the terminal degradation rate found with the DFOP model.")
        message("The representative half-life obtained from the DFOP model may be used")
      }
      t_rep <- t_DFOP2
    }
  }
  return(t_rep)
}