Skip to content

Commit

Permalink
fix: Uses tidyverse style coding
Browse files Browse the repository at this point in the history
* changes some argument names so they are self-describing and they use
  snake_case instead of . or camelcase
* it is no longer recommended to use return() at the end of a function
  if it is the last object, only if you are returning early
* removes some unused objects
  • Loading branch information
kellijohnson-NOAA committed Oct 6, 2024
1 parent ac638b9 commit 5335b8b
Show file tree
Hide file tree
Showing 6 changed files with 136 additions and 160 deletions.
220 changes: 121 additions & 99 deletions R/fimsfit.R
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,8 @@ is.FIMSFit <- function(x) inherits(x, "FIMSFit")
#' Check if an object is a list of FIMSFit objects
#' @param x List of fits returned from \code{\link{fit_tmb}}
#' @export
is.FIMSFits <- function(x){
if(!is.list(x)) {
is.FIMSFits <- function(x) {
if (!is.list(x)) {
cli::cli_warn(
message = c("x" = "{.par x} is not a list -- something went wrong.")
)
Expand All @@ -78,15 +78,18 @@ is.FIMSFits <- function(x){
#' @return Summary printed to console
#' @method print FIMSFit
#' @export
print.FIMSFit <- function(fit, ...){
rt <- as.numeric(fit@timing[["time_total"]], units='secs')
ru <- 'seconds'
if(rt>60*60*24) {
rt <- rt/(60*60*24); ru <- 'days'
} else if(rt>60*60) {
rt <- rt/(60*60); ru <- 'hours'
} else if(rt>60){
rt <- rt/60; ru <- 'minutes'
print.FIMSFit <- function(fit, ...) {
rt <- as.numeric(fit@timing[["time_total"]], units = "secs")
ru <- "seconds"
if (rt > 60 * 60 * 24) {
rt <- rt / (60 * 60 * 24)
ru <- "days"
} else if (rt > 60 * 60) {
rt <- rt / (60 * 60)
ru <- "hours"
} else if (rt > 60) {
rt <- rt / 60
ru <- "minutes"
}

number_of_parameters <- paste(
Expand Down Expand Up @@ -208,7 +211,7 @@ create_FIMSFit <- function(
estimates <- tibble::tibble(
as.data.frame(std)
) |>
dplyr::rename(value = Estimate, se = "Std. Error") |>
dplyr::rename(value = "Estimate", se = "Std. Error") |>
dplyr::mutate(
name = dimnames(std)[[1]],
.before = "value"
Expand All @@ -233,113 +236,132 @@ create_FIMSFit <- function(
timing = timing,
version = version
)
fit
}

#' Fit a FIMS model (BETA)
#' @param input Input list as returned by `prepare_input()`.
#' @param getsd Calculate and return sdreport?
#' @param loopnum A positive integer specifying the number of iterations of
#' the optimizer that will be performed to improve the gradient.
#' @param do.fit Optimize (TRUE, default) or (FALSE) build and return
#' @param get_sd Calculate and return sdreport?
#' @param number_of_loops A positive integer specifying the number of
#' iterations of the optimizer that will be performed to improve the
#' gradient.
#' @param optimize Optimize (TRUE, default) or (FALSE) build and return
#' a list containing the obj and report slot.
#' @param newtonsteps The number of Newton steps using the inverse
#' Hessian to do after optimization.
#' @param number_of_newton_steps The number of Newton steps using the inverse
#' Hessian to do after optimization. Not yet implemented.
#' @param control A list of optimizer settings passed to [stats::nlminb()].
#' @param filename Character string giving a file name to save
#' the fitted object as an RDS object. Defaults to 'fit.RDS',
#' and a value of NULL indicates not to save it. If specified,
#' it must end in .RDS. The file is written to folder given by
#' \code{input$path}.
#' \code{input$path}. Not yet implemented.
#' @return
#' An object of class `FIMSFit` is returned, where the structure is the same
#' regardless if `do.fit = TRUE` or not. More information is available in the
#' `fits` slot if `getsd = TRUE`, otherwise uncertainty information will not
#' regardless if `optimize = TRUE` or not. More information is available in the
#' `fits` slot if `get_sd = TRUE`, otherwise uncertainty information will not
#' be available.
#' @details This function is a beta version still and subject to change
#' without warning.
#' @export
fit_fims <- function(input, getsd=TRUE, loopnum=3, do.fit=TRUE, newtonsteps=0,
control=NULL,
filename=NULL){
if(!is.null(input$random)) {
cli::cli_abort("Random effects declared but not implemented yet.")
}
if(newtonsteps>0) {
cli::cli_abort("Newton steps not implemented yet.")
}
if (loopnum < 0) {
cli::cli_abort("loopnum ({.par {loopnum}}) must be >= 0.")
}
obj <- MakeADFun(data=list(), parameters=input$parameters,
map=input$map, random=input$random,
DLL='FIMS', silent=TRUE)
if(!do.fit) {
initial_fit <- create_FIMSFit(
input = input,
obj = obj,
timing = c("time_total" = as.difftime(0, units = "secs"))
fit_fims <- function(input,
get_sd = TRUE,
number_of_loops = 3,
optimize = TRUE,
number_of_newton_steps = 0,
control = NULL,
filename = NULL) {
if (!is.null(input$random)) {
cli::cli_abort("Random effects declared but not implemented yet.")
}
if (number_of_newton_steps > 0) {
cli::cli_abort("Newton steps not implemented yet.")
}
if (number_of_loops < 0) {
cli::cli_abort("number_of_loops ({.par {number_of_loops}}) must be >= 0.")
}
obj <- MakeADFun(
data = list(),
parameters = input$parameters,
map = input$map,
random = input$random,
DLL = "FIMS",
silent = TRUE
)
return(initial_fit)
}
# to do: max this update elements that are not supplied by default
if(is.null(control))
control <- list(eval.max=10000, iter.max=10000, trace=0)
if(!is_fims_verbose()) control$trace <- 0
## optimize and compare
cli::cli_inform(c("v" = "Starting optimization ..."))
t0 <- Sys.time()
opt0 <- opt <-
with(obj, nlminb(start = par, objective = fn, gradient = gr, control=control))
maxgrad0 <- maxgrad <- max(abs(obj$gr(opt$par)))
if(loopnum>0){
cli::cli_inform(c(
"i" = "Restarting optimizer {loopnum} times silently to improve gradient."
))
for(ii in 2:loopnum){
if (!optimize) {
initial_fit <- create_FIMSFit(
input = input,
obj = obj,
timing = c("time_total" = as.difftime(0, units = "secs"))
)
return(initial_fit)
}
# to do: max this update elements that are not supplied by default
if (is.null(control)) {
control <- list(eval.max = 10000, iter.max = 10000, trace = 0)
}
if (!is_fims_verbose()) {
control$trace <- 0
}
## optimize and compare
cli::cli_inform(c("v" = "Starting optimization ..."))
t0 <- Sys.time()
opt <- with(
obj,
nlminb(start = par, objective = fn, gradient = gr, control = control)
)
maxgrad0 <- maxgrad <- max(abs(obj$gr(opt$par)))
if (number_of_loops > 0) {
cli::cli_inform(c(
"i" = "Restarting optimizer {number_of_loops} times silently to improve gradient."
))
for (ii in 2:number_of_loops) {
# TODO: Shouldn't this be dictated by verbose?
control$trace <- 0
opt <- with(obj, nlminb(start = opt$par, objective = fn,
gradient = gr, control=control))
maxgrad <- max(abs(obj$gr(opt$par)))
opt <- with(
obj,
nlminb(start = opt$par, objective = fn, gradient = gr, control = control)
)
maxgrad <- max(abs(obj$gr(opt$par)))
}
div_digit <- cli::cli_div(theme = list(.val = list(digits = 5)))
cli::cli_inform(c(
"i" = "Maximum gradient went from {.val {maxgrad0}} to
{.val {maxgrad}} after {number_of_loops} steps."
))
cli::cli_end(div_digit)
}
div_digit <- cli::cli_div(theme = list(.val = list(digits = 5)))
cli::cli_inform(c(
"i" = "Maximum gradient went from {.val {maxgrad0}} to
{.val {maxgrad}} after {loopnum} steps."
))
cli::cli_end(div_digit)
}
time_optimization <- Sys.time() - t0
cli::cli_inform(c("v" = "Finished optimization"))
time_optimization <- Sys.time() - t0
cli::cli_inform(c("v" = "Finished optimization"))

sdrep <- std <- NULL
time_sdreport <- NA
if(getsd){
t2 <- Sys.time()
sd_report <- TMB::sdreport(obj)
cli::cli_inform(c("v" = "Finished sdreport"))
time_sdreport <- Sys.time() - t2
} else {
sd_report <- list()
time_sdreport <- as.difftime(0, units = "secs")
}
time_sdreport <- NA
if (get_sd) {
t2 <- Sys.time()
sd_report <- TMB::sdreport(obj)
cli::cli_inform(c("v" = "Finished sdreport"))
time_sdreport <- Sys.time() - t2
} else {
sd_report <- list()
time_sdreport <- as.difftime(0, units = "secs")
}

timing <- c(time_optimization=time_optimization,
time_sdreport=time_sdreport,
time_total=Sys.time() - t0)
fit <- create_FIMSFit(
input = input,
obj = obj,
opt = opt,
sd_report = sd_report,
timing = timing
)
print(fit)
if(!is.null(filename)) {
cli::cli_warn(c(
"i" = "Saving output to file is not yet implemented."
))
# saveRDS(fit, file=file.path(input[["path"]], filename))
}
return(fit)
timing <- c(
time_optimization = time_optimization,
time_sdreport = time_sdreport,
time_total = Sys.time() - t0
)
fit <- create_FIMSFit(
input = input,
obj = obj,
opt = opt,
sd_report = sd_report,
timing = timing
)
print(fit)
if (!is.null(filename)) {
cli::cli_warn(c(
"i" = "Saving output to file is not yet implemented."
))
# saveRDS(fit, file=file.path(input[["path"]], filename))
}
return(fit)
}
27 changes: 14 additions & 13 deletions man/fit_fims.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

14 changes: 0 additions & 14 deletions man/is.fimsfit.Rd

This file was deleted.

14 changes: 0 additions & 14 deletions man/is.fimsfits.Rd

This file was deleted.

19 changes: 0 additions & 19 deletions man/print.fimsfit.Rd

This file was deleted.

2 changes: 1 addition & 1 deletion vignettes/fims-demo.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -357,7 +357,7 @@ sucess <- CreateTMBModel()
parameters <- list(p = get_fixed())
input <- list(data = fims_frame, parameters = parameters, version = "FIMS demo")
# Run the model without optimization to help ensure a viable model
test_fit <- fit_fims(input, do.fit = FALSE)
test_fit <- fit_fims(input, optimize = FALSE)
# Run the model with optimization
fit <- fit_fims(input)
```
Expand Down

0 comments on commit 5335b8b

Please sign in to comment.