## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----runtimer, include=FALSE--------------------------------------------------
runstart <- lubridate::now()

## -----------------------------------------------------------------------------
# show all models in rTPC
rTPC::get_model_names()

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# eubank_1973 <- function(temp, topt, a, b) {
#   est <- a / ((temp - topt)^2 + b)
#   return(est)
# }

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# sharpeschoolhigh_1981 <- function(temp, r_tref, e, eh, th, tref) {
#   tref <- 273.15 + tref
#   k <- 8.62e-05
#   boltzmann.term <- r_tref * exp(e / k * (1 / tref - 1 / (temp + 273.15)))
#   inactivation.term <- 1 /
#     (1 + exp(eh / k * (1 / (th + 273.15) - 1 / (temp + 273.15))))
#   return(boltzmann.term * inactivation.term)
# }

## -----------------------------------------------------------------------------
# A slightly compressed example of how d is generated
library(rTPC)

subs <- subset(chlorella_tpc, chlorella_tpc$curve_id == 1)
d <- data.frame(x = subs$temp, y = subs$rate, stringsAsFactors = FALSE)
d <- d[order(d$x), ]
d

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# eubank_1973.starting_vals <- function(d) {
#   # starting values go here
#   return(c(topt = topt, a = a, b = b))
# }

## -----------------------------------------------------------------------------
eubank_1973.starting_vals <- function(d) {
  rmax = max(d$y, na.rm = TRUE) # Find max trait value
  topt = mean(d$x[d$y == rmax]) # Find T of rmax
  a = 300
  b = 50
  return(c(topt = topt, a = a, b = b))
}

## -----------------------------------------------------------------------------
eubank_1973.lower_lims <- function(d) {
  topt = 0
  a = 0
  b = 0
  return(c(topt = topt, a = a, b = b))
}

eubank_1973.upper_lims <- function(d) {
  topt = 150
  a = Inf
  b = Inf
  return(c(topt = topt, a = a, b = b))
}

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# # do not run the test on CRAN as they take too long
# testthat::skip_on_cran()
# 
# # method: fit model and get predictions. Check these against others.
# 
# # load in ggplot
# library(ggplot2)
# 
# # subset for the first TPC curve
# data('chlorella_tpc')
# d <- subset(chlorella_tpc, curve_id == 1)
# 
# # get start values and fit model
# start_vals <- get_start_vals(d$temp, d$rate, model_name = 'eubank_1973')
# 
# # fit model
# mod <- nls.multstart::nls_multstart(
#   rate ~ eubank_1973(temp = temp, tops, a, b),
#   data = d,
#   iter = c(3, 3, 3),
#   start_lower = start_vals - 10,
#   start_upper = start_vals + 10,
#   lower = get_lower_lims(d$temp, d$rate, model_name = 'eubank_1973'),
#   upper = get_upper_lims(d$temp, d$rate, model_name = 'eubank_1973'),
#   supp_errors = 'Y',
#   convergence_count = FALSE
# )
# 
# # get predictions
# preds <- broom::augment(mod)
# # dput(round(preds$.fitted, 1))
# 
# # plot
# ggplot(preds) +
#   geom_point(aes(temp, rate)) +
#   geom_line(aes(temp, .fitted)) +
#   theme_bw()
# 
# # run test
# testthat::test_that("eubank_1973 function works", {
#   testthat::expect_equal(
#     round(preds$.fitted, 1),
#     c(0.2, 0.2, 0.3, 0.4, 0.6, 0.9, 1.3, 1.6, 1.4, 1, 0.7, 0.4)
#   )
# })

## ----eval=FALSE, echo=TRUE----------------------------------------------------
# get_start_vals <- function(x, y, model_name) {
#   mod_names <- get_model_names(returnall = TRUE)
#   model_name <- tryCatch(
#     rlang::arg_match(model_name, mod_names),
#     error = function(e) {
#       cli::cli_abort(
#         c(
#           "x" = "Supplied {.arg model_name} ({.val {model_name}}) is not an available model in rTPC.",
#           "!" = "Please check the spelling of {.arg model_name}.",
#           " " = "(run {.fn rTPC::get_model_names} to see all valid names.)",
#           ""
#         ),
#         call = rlang::caller_env(n = 4)
#       )
#     }
#   )
# 
#   # make data frame
#   d <- data.frame(x, y, stringsAsFactors = FALSE)
#   d <- d[order(d$x), ]
# 
#   start_vals <- tryCatch(
#     do.call(paste0(model_name, ".starting_vals"), list(d = d)),
#     error = function(e) {
#       NULL
#     }
#   )
# 
#   return(start_vals)
# }

## -----------------------------------------------------------------------------
model_name <- "eubank_1976"
paste0(model_name, ".starting_vals")

## ----tot_time, include=FALSE--------------------------------------------------
tot_time <- lubridate::as.duration(lubridate::now() - runstart)