R/pr_fit.R

In khufkens/phenor: An R phenology Modelling Framework

#' Fit a single phenology model
#'
#' General purpose model fit function, called by higher order
#' pr_fit_comaprison() and pr_cross_validate() functions. In itself a wrapper
#' around the pr_fit_parameter() function.
#'
#' This function together with pr_predict() is in many ways
#' the core of the package.
#'
#' @param model the model name to be used in optimizing the model
#' @param data dataset generated by the format_phenocam() or
#' format_modis() routines, or adhering to the general model optimization
#' input format.
#' @param method optimization method to use (default = GenSA)
#'    - GenSA :  Generalized Simulated Annealing algorithm
#'    - genoud : GENetic Optimization Using Derivatives
#'    - BayesianTools: various bayesian based optimization tools
#' @param par_ranges a vector of starting parameter values (function specific)
#' defaults to the parameter ranges as provided with the current models
#' and set forth by Basler (2016)
#' @param control list of control parameters to be passed to the optimizer
#' @param plot TRUE / FALSE, plot model fit
#' @param ... additional control parameters to be passed
#' @keywords phenology, model, validation
#' @export
#' @examples
#'
#' \dontrun{
#' pr_fit(model,
#'        par_ranges = "parameter_ranges.csv")
#' }
#This is a test

pr_fit <- function(
  model = "TT",
  data = phenor::phenocam_DB,
  method = "GenSA",
  control = list(max.call = 2000),
  par_ranges = system.file(
    "extdata",
    "parameter_ranges.csv",
    package = "phenor",
    mustWork = TRUE
  ),
  plot = TRUE,
  ...
  ){

  # convert to a flat format for speed
  data <- pr_flatten(data)

  # read parameter ranges
  d <- pr_parameters(model = model,
                     par_ranges = par_ranges)

  # optimize paramters
  optim_par = pr_fit_parameters(
    par = NULL,
    data = data,
    model = model,
    method = method,
    lower = as.numeric(d[1,]),
    upper = as.numeric(d[2,]),
    control = control,
    ...
  )

  # estimate model output using optimized
  # parameters
  out <- pr_predict(
    data = data,
    model = model,
    par = optim_par$par
  )

  # basic statistics
  RMSE_NULL <- sqrt(mean((data$transition_dates - null(data)) ^ 2, na.rm = T))
  RMSE <- rmse(par = optim_par$par, data = data, model = model)
  Ac <- AICc(measured = data$transition_dates,
            predicted = out,
            k = length(optim_par$par))

  output <- list(
    "model" = model,
    "par" = optim_par$par,
    "measured" = data$transition_dates,
    "predicted" = out,
    "rmse" = RMSE,
    "rmse_null" = RMSE_NULL,
    "aic" = Ac,
    "opt_out" = optim_par$opt_out
  )

  # assign class
  class(output) <- "pr_fit"

  # return optimized parameters and stats
  return(output)
}