From a11faadb22713de700f181e3464ec780b2883a9a Mon Sep 17 00:00:00 2001 From: ChristineStawitz-NOAA Date: Tue, 28 Nov 2023 18:38:06 +0000 Subject: [PATCH] style and docs: run devtools::document() and styler::style_pkg() --- R/fimsframe.R | 8 +- tests/concurrent/fims_concurrent_mpi.R | 387 ++++++++++---------- tests/concurrent/fims_concurrent_processR.R | 123 +++---- tests/concurrent/fims_concurrent_snowfall.R | 360 +++++++++--------- tests/testthat/test-data-object.R | 2 - tests/testthat/test-fims-estimation.R | 32 +- tests/testthat/test-fims-rcpp.R | 1 - tests/testthat/test-fleet-interface.R | 1 - tests/testthat/test-maturity-interface.R | 1 - tests/testthat/test-recruitment-interface.R | 1 - tests/testthat/test-selectivity-interface.R | 1 - vignettes/fims-demo.Rmd | 14 +- 12 files changed, 461 insertions(+), 470 deletions(-) diff --git a/R/fimsframe.R b/R/fimsframe.R index 85535ba6..0efd4c0a 100644 --- a/R/fimsframe.R +++ b/R/fimsframe.R @@ -16,7 +16,7 @@ setClass( data = "data.frame", # can use c( ) or list here. fleets = "numeric", nyrs = "integer", - start_year = "integer", + start_year = "integer", end_year = "integer" ) ) @@ -260,7 +260,7 @@ setValidity( if (!"dateend" %in% colnames(object@data)) { errors <- c(errors, "data must contain 'uncertainty'") } - + # TODO: Add checks for other slots @@ -327,7 +327,7 @@ FIMSFrame <- function(data) { data = data, fleets = fleets, nyrs = nyrs, - start_year = start_year, + start_year = start_year, end_year = end_year ) return(out) @@ -357,7 +357,7 @@ FIMSFrameAge <- function(data) { data = data, fleets = fleets, nyrs = nyrs, - start_year = start_year, + start_year = start_year, end_year = end_year, ages = ages, nages = nages, diff --git a/tests/concurrent/fims_concurrent_mpi.R b/tests/concurrent/fims_concurrent_mpi.R index 2ae702ed..bd15bf29 100644 --- a/tests/concurrent/fims_concurrent_mpi.R +++ b/tests/concurrent/fims_concurrent_mpi.R @@ -3,21 +3,21 @@ library(minimizR) # Load the R MPI package if it is not already loaded. if (!is.loaded("mpi_initialize")) { - library("Rmpi") + library("Rmpi") } # # In case R exits unexpectedly, have it automatically clean up # resources taken up by Rmpi (slaves, memory, etc...) .Last <- function() { - if (is.loaded("mpi_initialize")) { - if (mpi.comm.size(1) > 0) { - print("Please use mpi.close.Rslaves() to close slaves.") - mpi.close.Rslaves() - } - print("Please use mpi.quit() to quit R") - .Call("mpi_finalize") + if (is.loaded("mpi_initialize")) { + if (mpi.comm.size(1) > 0) { + print("Please use mpi.close.Rslaves() to close slaves.") + mpi.close.Rslaves() } + print("Please use mpi.quit() to quit R") + .Call("mpi_finalize") + } } @@ -27,14 +27,14 @@ working_dir <- getwd() maindir <- tempdir() model_input <- ASSAMC::save_initial_input() FIMS_C0_estimation <- ASSAMC::save_initial_input( -base_case = TRUE, -input_list = model_input, -maindir = maindir, -om_sim_num = 1, -keep_sim_num = 1, -figure_number = 1, -seed_num = 9924, -case_name = "FIMS_C0_estimation" + base_case = TRUE, + input_list = model_input, + maindir = maindir, + om_sim_num = 1, + keep_sim_num = 1, + figure_number = 1, + seed_num = 9924, + case_name = "FIMS_C0_estimation" ) ASSAMC::run_om(input_list = FIMS_C0_estimation) @@ -45,166 +45,166 @@ setwd(working_dir) -NUMBER_OF_MODEL_RUNS<-100 +NUMBER_OF_MODEL_RUNS <- 100 # Fix and change R0 randomly to get different output for each # model run. -init_fims<-function(i){ - fims <- Rcpp::Module("fims", PACKAGE = "FIMS") - - # Recruitment - recruitment <- new(fims$BevertonHoltRecruitment) - # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is taken - # before the likelihood calculation - recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) - recruitment$log_rzero$value <- log(om_input$R0 + runif(1,min=0, max=1000)) - recruitment$log_rzero$is_random_effect <- FALSE - recruitment$log_rzero$estimated <- FALSE - recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) - recruitment$logit_steep$is_random_effect <- FALSE - recruitment$logit_steep$estimated <- FALSE - recruitment$estimate_deviations <- TRUE - # recruit deviations should enter the model in normal space. - # The log is taken in the likelihood calculations - # alternative setting: recruitment$deviations <- rep(1, length(om_input$logR.resid)) - recruitment$deviations <- exp(om_input$logR.resid) - - # Data - catch <- em_input$L.obs$fleet1 - fishing_fleet_index <- new(fims$Index, length(catch)) - fishing_fleet_index$index_data <- catch - fishing_fleet_age_comp <- new(fims$AgeComp, length(catch), om_input$nages) - fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1 - - survey_index <- em_input$surveyB.obs$survey1 - survey_fleet_index <- new(fims$Index, length(survey_index)) - survey_fleet_index$index_data <- survey_index - survey_fleet_age_comp <- new(fims$AgeComp, length(survey_index), om_input$nages) - survey_fleet_age_comp$age_comp_data <- c(t(em_input$survey.age.obs$survey1)) * em_input$n.survey$survey1 - - # Growth - ewaa_growth <- new(fims$EWAAgrowth) - ewaa_growth$ages <- om_input$ages - ewaa_growth$weights <- om_input$W.mt - - # Maturity - maturity <- new(fims$LogisticMaturity) - maturity$median$value <- om_input$A50.mat - maturity$median$is_random_effect <- FALSE - maturity$median$estimated <- FALSE - maturity$slope$value <- om_input$slope - maturity$slope$is_random_effect <- FALSE - maturity$slope$estimated <- FALSE - - # Fleet - # Create the fishing fleet - fishing_fleet_selectivity <- new(fims$LogisticSelectivity) - fishing_fleet_selectivity$median$value <- om_input$sel_fleet$fleet1$A50.sel1 - fishing_fleet_selectivity$median$is_random_effect <- FALSE - fishing_fleet_selectivity$median$estimated <- TRUE - fishing_fleet_selectivity$slope$value <- om_input$sel_fleet$fleet1$slope.sel1 - fishing_fleet_selectivity$slope$is_random_effect <- FALSE - fishing_fleet_selectivity$slope$estimated <- TRUE - - fishing_fleet <- new(fims$Fleet) - fishing_fleet$nages <- om_input$nages - fishing_fleet$nyears <- om_input$nyr - fishing_fleet$log_Fmort <- log(om_output$f) - fishing_fleet$estimate_F <- TRUE - fishing_fleet$random_F <- FALSE - fishing_fleet$log_q <- log(1.0) - fishing_fleet$estimate_q <- FALSE - fishing_fleet$random_q <- FALSE - fishing_fleet$log_obs_error$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1))) - fishing_fleet$log_obs_error$estimated <- FALSE - # Need get_id() for setting up observed agecomp and index data? - fishing_fleet$SetAgeCompLikelihood(1) - fishing_fleet$SetIndexLikelihood(1) - fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id()) - fishing_fleet$SetObservedIndexData(fishing_fleet_index$get_id()) - fishing_fleet$SetObservedAgeCompData(fishing_fleet_age_comp$get_id()) - - # Create the survey fleet - survey_fleet_selectivity <- new(fims$LogisticSelectivity) - survey_fleet_selectivity$median$value <- om_input$sel_survey$survey1$A50.sel1 - survey_fleet_selectivity$median$is_random_effect <- FALSE - survey_fleet_selectivity$median$estimated <- TRUE - survey_fleet_selectivity$slope$value <- om_input$sel_survey$survey1$slope.sel1 - survey_fleet_selectivity$slope$is_random_effect <- FALSE - survey_fleet_selectivity$slope$estimated <- TRUE - - survey_fleet <- new(fims$Fleet) - survey_fleet$is_survey <- TRUE - survey_fleet$nages <- om_input$nages - survey_fleet$nyears <- om_input$nyr - # survey_fleet$log_Fmort <- rep(log(0.0000000000000000000000000001), om_input$nyr) #-Inf? - survey_fleet$estimate_F <- FALSE - survey_fleet$random_F <- FALSE - survey_fleet$log_q <- log(om_output$survey_q$survey1) - survey_fleet$estimate_q <- TRUE - survey_fleet$random_q <- FALSE - survey_fleet$log_obs_error$value <- log(sqrt(log(em_input$cv.survey$survey1^2 + 1))) - survey_fleet$log_obs_error$estimated <- FALSE - survey_fleet$SetAgeCompLikelihood(1) - survey_fleet$SetIndexLikelihood(1) - survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id()) - survey_fleet$SetObservedIndexData(survey_fleet_index$get_id()) - survey_fleet$SetObservedAgeCompData(survey_fleet_age_comp$get_id()) - - # Population - population <- new(fims$Population) - # is it a problem these are not Parameters in the Population interface? - # the Parameter class (from rcpp/rcpp_objects/rcpp_interface_base) cannot handle vectors, - # do we need a ParameterVector class? - population$log_M <- rep(log(om_input$M.age[1]), om_input$nyr * om_input$nages) - population$estimate_M <- FALSE - population$log_init_naa <- log(om_output$N.age[1, ]) - population$estimate_init_naa <- TRUE - population$nages <- om_input$nages - population$ages <- om_input$ages - population$nfleets <- sum(om_input$fleet_num, om_input$survey_num) - population$nseasons <- 1 - population$nyears <- om_input$nyr - population$prop_female <- om_input$proportion.female[1] - population$SetMaturity(maturity$get_id()) - population$SetGrowth(ewaa_growth$get_id()) - population$SetRecruitment(recruitment$get_id()) - - ## Set-up TMB - fims$CreateTMBModel() - parameters <- list(p = fims$get_fixed()) - obj<-TMB::MakeADFun(data = list(), parameters, DLL = "FIMS") - - fims$clear() - return(obj) +init_fims <- function(i) { + fims <- Rcpp::Module("fims", PACKAGE = "FIMS") + + # Recruitment + recruitment <- new(fims$BevertonHoltRecruitment) + # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is taken + # before the likelihood calculation + recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) + recruitment$log_rzero$value <- log(om_input$R0 + runif(1, min = 0, max = 1000)) + recruitment$log_rzero$is_random_effect <- FALSE + recruitment$log_rzero$estimated <- FALSE + recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) + recruitment$logit_steep$is_random_effect <- FALSE + recruitment$logit_steep$estimated <- FALSE + recruitment$estimate_deviations <- TRUE + # recruit deviations should enter the model in normal space. + # The log is taken in the likelihood calculations + # alternative setting: recruitment$deviations <- rep(1, length(om_input$logR.resid)) + recruitment$deviations <- exp(om_input$logR.resid) + + # Data + catch <- em_input$L.obs$fleet1 + fishing_fleet_index <- new(fims$Index, length(catch)) + fishing_fleet_index$index_data <- catch + fishing_fleet_age_comp <- new(fims$AgeComp, length(catch), om_input$nages) + fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1 + + survey_index <- em_input$surveyB.obs$survey1 + survey_fleet_index <- new(fims$Index, length(survey_index)) + survey_fleet_index$index_data <- survey_index + survey_fleet_age_comp <- new(fims$AgeComp, length(survey_index), om_input$nages) + survey_fleet_age_comp$age_comp_data <- c(t(em_input$survey.age.obs$survey1)) * em_input$n.survey$survey1 + + # Growth + ewaa_growth <- new(fims$EWAAgrowth) + ewaa_growth$ages <- om_input$ages + ewaa_growth$weights <- om_input$W.mt + + # Maturity + maturity <- new(fims$LogisticMaturity) + maturity$median$value <- om_input$A50.mat + maturity$median$is_random_effect <- FALSE + maturity$median$estimated <- FALSE + maturity$slope$value <- om_input$slope + maturity$slope$is_random_effect <- FALSE + maturity$slope$estimated <- FALSE + + # Fleet + # Create the fishing fleet + fishing_fleet_selectivity <- new(fims$LogisticSelectivity) + fishing_fleet_selectivity$median$value <- om_input$sel_fleet$fleet1$A50.sel1 + fishing_fleet_selectivity$median$is_random_effect <- FALSE + fishing_fleet_selectivity$median$estimated <- TRUE + fishing_fleet_selectivity$slope$value <- om_input$sel_fleet$fleet1$slope.sel1 + fishing_fleet_selectivity$slope$is_random_effect <- FALSE + fishing_fleet_selectivity$slope$estimated <- TRUE + + fishing_fleet <- new(fims$Fleet) + fishing_fleet$nages <- om_input$nages + fishing_fleet$nyears <- om_input$nyr + fishing_fleet$log_Fmort <- log(om_output$f) + fishing_fleet$estimate_F <- TRUE + fishing_fleet$random_F <- FALSE + fishing_fleet$log_q <- log(1.0) + fishing_fleet$estimate_q <- FALSE + fishing_fleet$random_q <- FALSE + fishing_fleet$log_obs_error$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1))) + fishing_fleet$log_obs_error$estimated <- FALSE + # Need get_id() for setting up observed agecomp and index data? + fishing_fleet$SetAgeCompLikelihood(1) + fishing_fleet$SetIndexLikelihood(1) + fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id()) + fishing_fleet$SetObservedIndexData(fishing_fleet_index$get_id()) + fishing_fleet$SetObservedAgeCompData(fishing_fleet_age_comp$get_id()) + + # Create the survey fleet + survey_fleet_selectivity <- new(fims$LogisticSelectivity) + survey_fleet_selectivity$median$value <- om_input$sel_survey$survey1$A50.sel1 + survey_fleet_selectivity$median$is_random_effect <- FALSE + survey_fleet_selectivity$median$estimated <- TRUE + survey_fleet_selectivity$slope$value <- om_input$sel_survey$survey1$slope.sel1 + survey_fleet_selectivity$slope$is_random_effect <- FALSE + survey_fleet_selectivity$slope$estimated <- TRUE + + survey_fleet <- new(fims$Fleet) + survey_fleet$is_survey <- TRUE + survey_fleet$nages <- om_input$nages + survey_fleet$nyears <- om_input$nyr + # survey_fleet$log_Fmort <- rep(log(0.0000000000000000000000000001), om_input$nyr) #-Inf? + survey_fleet$estimate_F <- FALSE + survey_fleet$random_F <- FALSE + survey_fleet$log_q <- log(om_output$survey_q$survey1) + survey_fleet$estimate_q <- TRUE + survey_fleet$random_q <- FALSE + survey_fleet$log_obs_error$value <- log(sqrt(log(em_input$cv.survey$survey1^2 + 1))) + survey_fleet$log_obs_error$estimated <- FALSE + survey_fleet$SetAgeCompLikelihood(1) + survey_fleet$SetIndexLikelihood(1) + survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id()) + survey_fleet$SetObservedIndexData(survey_fleet_index$get_id()) + survey_fleet$SetObservedAgeCompData(survey_fleet_age_comp$get_id()) + + # Population + population <- new(fims$Population) + # is it a problem these are not Parameters in the Population interface? + # the Parameter class (from rcpp/rcpp_objects/rcpp_interface_base) cannot handle vectors, + # do we need a ParameterVector class? + population$log_M <- rep(log(om_input$M.age[1]), om_input$nyr * om_input$nages) + population$estimate_M <- FALSE + population$log_init_naa <- log(om_output$N.age[1, ]) + population$estimate_init_naa <- TRUE + population$nages <- om_input$nages + population$ages <- om_input$ages + population$nfleets <- sum(om_input$fleet_num, om_input$survey_num) + population$nseasons <- 1 + population$nyears <- om_input$nyr + population$prop_female <- om_input$proportion.female[1] + population$SetMaturity(maturity$get_id()) + population$SetGrowth(ewaa_growth$get_id()) + population$SetRecruitment(recruitment$get_id()) + + ## Set-up TMB + fims$CreateTMBModel() + parameters <- list(p = fims$get_fixed()) + obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS") + + fims$clear() + return(obj) } -run_fims<-function(begin,end){ - library(FIMS) - library(minimizR) - minimizer<-1 - results<-list() - - id <- mpi.comm.rank() - index<-1 - for(i in begin[id]:end[id]){ - obj<-init_fims(i) - if(minimizer == 1){ - results[[index]] <- optim(obj$par, obj$fn, obj$gr, - method = "BFGS", - control = list(maxit = 1000000, reltol = 1e-15) - ) - }else{ - results[[index]]<-minimizR(obj$par, obj$fn, obj$gr, control = list(tolerance = 1e-8)) - } - - index = index+1 +run_fims <- function(begin, end) { + library(FIMS) + library(minimizR) + minimizer <- 1 + results <- list() + + id <- mpi.comm.rank() + index <- 1 + for (i in begin[id]:end[id]) { + obj <- init_fims(i) + if (minimizer == 1) { + results[[index]] <- optim(obj$par, obj$fn, obj$gr, + method = "BFGS", + control = list(maxit = 1000000, reltol = 1e-15) + ) + } else { + results[[index]] <- minimizR(obj$par, obj$fn, obj$gr, control = list(tolerance = 1e-8)) } - return(results) + + index <- index + 1 + } + return(results) } @@ -221,51 +221,50 @@ nsims <- NUMBER_OF_MODEL_RUNS begin <- rep(0, ns) end <- rep(0, ns) -#create scenario segments +# create scenario segments if (id == 0) { - segments <- nsims / ns - print(paste("segments ", segments)) - for (i in 1:ns) { - if (i < ns) { - begin[i] <- as.integer((i - 1) * segments+1) - end[i] <- as.integer(i * segments) - } else{ - begin[i] <- as.integer((i - 1) * segments+1) - end[i] <- nsims - } + segments <- nsims / ns + print(paste("segments ", segments)) + for (i in 1:ns) { + if (i < ns) { + begin[i] <- as.integer((i - 1) * segments + 1) + end[i] <- as.integer(i * segments) + } else { + begin[i] <- as.integer((i - 1) * segments + 1) + end[i] <- nsims } - print(begin) - print(end) + } + print(begin) + print(end) } mpi.spawn.Rslaves(nslaves = ns) -#capture minimizer output -completed<-list() +# capture minimizer output +completed <- list() -#pass the function to all slaves -mpi.bcast.Robj2slave( all = TRUE) -#mpi.bcast.Robj2slave(obj = run_fims) -#mpi.bcast.Robj2slave(obj = init_fims) +# pass the function to all slaves +mpi.bcast.Robj2slave(all = TRUE) +# mpi.bcast.Robj2slave(obj = run_fims) +# mpi.bcast.Robj2slave(obj = init_fims) -start<-Sys.time() -#execute all slaves and append the completed list -append(completed,mpi.remote.exec(run_fims, begin, end)) -end<-Sys.time() +start <- Sys.time() +# execute all slaves and append the completed list +append(completed, mpi.remote.exec(run_fims, begin, end)) +end <- Sys.time() print(completed) -runtime<- (end - start) -print(paste0(paste0(paste0(NUMBER_OF_MODEL_RUNS," model runs completed in "),runtime)," seconds.")) +runtime <- (end - start) +print(paste0(paste0(paste0(NUMBER_OF_MODEL_RUNS, " model runs completed in "), runtime), " seconds.")) -line=paste("Rmpi runtime ", runtime) -write(line,file="time.txt",append=TRUE) +line <- paste("Rmpi runtime ", runtime) +write(line, file = "time.txt", append = TRUE) # Tell all slaves to close down, and exit the program mpi.close.Rslaves(dellog = FALSE) mpi.quit() - diff --git a/tests/concurrent/fims_concurrent_processR.R b/tests/concurrent/fims_concurrent_processR.R index bfc9332d..9aceb6fb 100644 --- a/tests/concurrent/fims_concurrent_processR.R +++ b/tests/concurrent/fims_concurrent_processR.R @@ -3,7 +3,7 @@ library(minimizR) library(processR) -#get the Rcpp module +# get the Rcpp module p <- Rcpp::Module(module = "processR", PACKAGE = "processR") @@ -32,20 +32,20 @@ setwd(working_dir) -NUMBER_OF_MODEL_RUNS<-100 +NUMBER_OF_MODEL_RUNS <- 100 # Fix and change R0 randomly to get different output for each # model run. -init_fims<-function(i){ +init_fims <- function(i) { fims <- Rcpp::Module("fims", PACKAGE = "FIMS") - + # Recruitment recruitment <- new(fims$BevertonHoltRecruitment) # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is taken # before the likelihood calculation recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) - recruitment$log_rzero$value <- log(om_input$R0 + runif(1,min=0, max=1000)) + recruitment$log_rzero$value <- log(om_input$R0 + runif(1, min = 0, max = 1000)) recruitment$log_rzero$is_random_effect <- FALSE recruitment$log_rzero$estimated <- FALSE recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) @@ -56,25 +56,25 @@ init_fims<-function(i){ # The log is taken in the likelihood calculations # alternative setting: recruitment$deviations <- rep(1, length(om_input$logR.resid)) recruitment$deviations <- exp(om_input$logR.resid) - + # Data catch <- em_input$L.obs$fleet1 fishing_fleet_index <- new(fims$Index, length(catch)) fishing_fleet_index$index_data <- catch fishing_fleet_age_comp <- new(fims$AgeComp, length(catch), om_input$nages) fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1 - + survey_index <- em_input$surveyB.obs$survey1 survey_fleet_index <- new(fims$Index, length(survey_index)) survey_fleet_index$index_data <- survey_index survey_fleet_age_comp <- new(fims$AgeComp, length(survey_index), om_input$nages) survey_fleet_age_comp$age_comp_data <- c(t(em_input$survey.age.obs$survey1)) * em_input$n.survey$survey1 - + # Growth ewaa_growth <- new(fims$EWAAgrowth) ewaa_growth$ages <- om_input$ages ewaa_growth$weights <- om_input$W.mt - + # Maturity maturity <- new(fims$LogisticMaturity) maturity$median$value <- om_input$A50.mat @@ -83,7 +83,7 @@ init_fims<-function(i){ maturity$slope$value <- om_input$slope maturity$slope$is_random_effect <- FALSE maturity$slope$estimated <- FALSE - + # Fleet # Create the fishing fleet fishing_fleet_selectivity <- new(fims$LogisticSelectivity) @@ -93,7 +93,7 @@ init_fims<-function(i){ fishing_fleet_selectivity$slope$value <- om_input$sel_fleet$fleet1$slope.sel1 fishing_fleet_selectivity$slope$is_random_effect <- FALSE fishing_fleet_selectivity$slope$estimated <- TRUE - + fishing_fleet <- new(fims$Fleet) fishing_fleet$nages <- om_input$nages fishing_fleet$nyears <- om_input$nyr @@ -111,7 +111,7 @@ init_fims<-function(i){ fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id()) fishing_fleet$SetObservedIndexData(fishing_fleet_index$get_id()) fishing_fleet$SetObservedAgeCompData(fishing_fleet_age_comp$get_id()) - + # Create the survey fleet survey_fleet_selectivity <- new(fims$LogisticSelectivity) survey_fleet_selectivity$median$value <- om_input$sel_survey$survey1$A50.sel1 @@ -120,7 +120,7 @@ init_fims<-function(i){ survey_fleet_selectivity$slope$value <- om_input$sel_survey$survey1$slope.sel1 survey_fleet_selectivity$slope$is_random_effect <- FALSE survey_fleet_selectivity$slope$estimated <- TRUE - + survey_fleet <- new(fims$Fleet) survey_fleet$is_survey <- TRUE survey_fleet$nages <- om_input$nages @@ -138,7 +138,7 @@ init_fims<-function(i){ survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id()) survey_fleet$SetObservedIndexData(survey_fleet_index$get_id()) survey_fleet$SetObservedAgeCompData(survey_fleet_age_comp$get_id()) - + # Population population <- new(fims$Population) # is it a problem these are not Parameters in the Population interface? @@ -157,12 +157,12 @@ init_fims<-function(i){ population$SetMaturity(maturity$get_id()) population$SetGrowth(ewaa_growth$get_id()) population$SetRecruitment(recruitment$get_id()) - + ## Set-up TMB fims$CreateTMBModel() parameters <- list(p = fims$get_fixed()) - obj<-TMB::MakeADFun(data = list(), parameters, DLL = "FIMS") - + obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS") + fims$clear() return(obj) } @@ -171,28 +171,28 @@ init_fims<-function(i){ -run_fims<-function(){ - results<<-list() +run_fims <- function() { + results <<- list() library(FIMS) library(minimizR) - minimizer<-2 + minimizer <- 2 + - id <- processR.rank - index<-1 - for(i in begin[id]:end[id]){ + index <- 1 + for (i in begin[id]:end[id]) { print(paste(id, paste(i, "me be here"))) - obj<-init_fims(i) - if(minimizer == 1){ - results[[index]] <<- optim(obj$par, obj$fn, obj$gr, - method = "BFGS", - control = list(maxit = 1000000, reltol = 1e-15) + obj <- init_fims(i) + if (minimizer == 1) { + results[[index]] <<- optim(obj$par, obj$fn, obj$gr, + method = "BFGS", + control = list(maxit = 1000000, reltol = 1e-15) ) - }else{ - results[[index]]<<-minimizR(obj$par, obj$fn, obj$gr, control = list(tolerance = 1e-4, verbose = FALSE)) + } else { + results[[index]] <<- minimizR(obj$par, obj$fn, obj$gr, control = list(tolerance = 1e-4, verbose = FALSE)) } - #write(x = i, file = paste0(paste("id",processR.rank)),".txt") - index = index+1 + # write(x = i, file = paste0(paste("id",processR.rank)),".txt") + index <- index + 1 } # return(results) } @@ -203,24 +203,24 @@ run_fims<-function(){ -id <- 0#mpi.comm.rank(comm = 0) -ns <- processR::HardwareConcurrency() +id <- 0 # mpi.comm.rank(comm = 0) +ns <- processR::HardwareConcurrency() nsims <- NUMBER_OF_MODEL_RUNS begin <- rep(0, ns) end <- rep(0, ns) -#create scenario segments +# create scenario segments if (id == 0) { segments <- nsims / ns print(paste("segments ", segments)) for (i in 1:ns) { if (i < ns) { - begin[i] <- as.integer((i - 1) * segments+1) + begin[i] <- as.integer((i - 1) * segments + 1) end[i] <- as.integer(i * segments) - } else{ - begin[i] <- as.integer((i - 1) * segments+1) + } else { + begin[i] <- as.integer((i - 1) * segments + 1) end[i] <- nsims } } @@ -230,43 +230,42 @@ if (id == 0) { -start<-Sys.time() +start <- Sys.time() -#create a pool of child processes -pool<-list() +# create a pool of child processes +pool <- list() -for(i in 1:processR::HardwareConcurrency()){ - #create a new child process +for (i in 1:processR::HardwareConcurrency()) { + # create a new child process pool[[i]] <- new(p$Process) - - #start the process. pass the function, environment, and rank + + # start the process. pass the function, environment, and rank pool[[i]]$start(run_fims, environment(), i) } -#iterate of the children and capture information -for(i in 1:processR::HardwareConcurrency()){ - - #wait for the child to complete - pool[[i]]$wait() - - #get child out stream - #message<-pool[[i]]$get_message() - - #access the childs environment - env<-as.environment(pool[[i]]$get_environment()); - #show minimizer results +# iterate of the children and capture information +for (i in 1:processR::HardwareConcurrency()) { + # wait for the child to complete + pool[[i]]$wait() + + # get child out stream + # message<-pool[[i]]$get_message() + + # access the childs environment + env <- as.environment(pool[[i]]$get_environment()) + # show minimizer results # print(env[["results"]]) } -end_<-Sys.time() +end_ <- Sys.time() -runtime<- (end_ - start) -print(paste0(paste0(paste0(NUMBER_OF_MODEL_RUNS," model runs completed in "),runtime)," seconds.")) +runtime <- (end_ - start) +print(paste0(paste0(paste0(NUMBER_OF_MODEL_RUNS, " model runs completed in "), runtime), " seconds.")) print(begin) print(end) -line=paste("processR runtime ", runtime) -write(line,file="time.txt",append=TRUE) +line <- paste("processR runtime ", runtime) +write(line, file = "time.txt", append = TRUE) diff --git a/tests/concurrent/fims_concurrent_snowfall.R b/tests/concurrent/fims_concurrent_snowfall.R index 1074a69c..93d5cc10 100644 --- a/tests/concurrent/fims_concurrent_snowfall.R +++ b/tests/concurrent/fims_concurrent_snowfall.R @@ -9,14 +9,14 @@ working_dir <- getwd() maindir <- tempdir() model_input <- ASSAMC::save_initial_input() FIMS_C0_estimation <- ASSAMC::save_initial_input( -base_case = TRUE, -input_list = model_input, -maindir = maindir, -om_sim_num = 1, -keep_sim_num = 1, -figure_number = 1, -seed_num = 9924, -case_name = "FIMS_C0_estimation" + base_case = TRUE, + input_list = model_input, + maindir = maindir, + om_sim_num = 1, + keep_sim_num = 1, + figure_number = 1, + seed_num = 9924, + case_name = "FIMS_C0_estimation" ) ASSAMC::run_om(input_list = FIMS_C0_estimation) @@ -27,166 +27,166 @@ setwd(working_dir) -NUMBER_OF_MODEL_RUNS<-100 +NUMBER_OF_MODEL_RUNS <- 100 # Fix and change R0 randomly to get different output for each # model run. -init_fims<-function(i){ - fims <- Rcpp::Module("fims", PACKAGE = "FIMS") - - # Recruitment - recruitment <- new(fims$BevertonHoltRecruitment) - # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is taken - # before the likelihood calculation - recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) - recruitment$log_rzero$value <- log(om_input$R0)# + runif(1,min=0, max=1000)) - recruitment$log_rzero$is_random_effect <- FALSE - recruitment$log_rzero$estimated <- TRUE - recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) - recruitment$logit_steep$is_random_effect <- FALSE - recruitment$logit_steep$estimated <- FALSE - recruitment$estimate_deviations <- TRUE - # recruit deviations should enter the model in normal space. - # The log is taken in the likelihood calculations - # alternative setting: recruitment$deviations <- rep(1, length(om_input$logR.resid)) - recruitment$deviations <- exp(om_input$logR.resid) - - # Data - catch <- em_input$L.obs$fleet1 - fishing_fleet_index <- new(fims$Index, length(catch)) - fishing_fleet_index$index_data <- catch - fishing_fleet_age_comp <- new(fims$AgeComp, length(catch), om_input$nages) - fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1 - - survey_index <- em_input$surveyB.obs$survey1 - survey_fleet_index <- new(fims$Index, length(survey_index)) - survey_fleet_index$index_data <- survey_index - survey_fleet_age_comp <- new(fims$AgeComp, length(survey_index), om_input$nages) - survey_fleet_age_comp$age_comp_data <- c(t(em_input$survey.age.obs$survey1)) * em_input$n.survey$survey1 - - # Growth - ewaa_growth <- new(fims$EWAAgrowth) - ewaa_growth$ages <- om_input$ages - ewaa_growth$weights <- om_input$W.mt - - # Maturity - maturity <- new(fims$LogisticMaturity) - maturity$median$value <- om_input$A50.mat - maturity$median$is_random_effect <- FALSE - maturity$median$estimated <- FALSE - maturity$slope$value <- om_input$slope - maturity$slope$is_random_effect <- FALSE - maturity$slope$estimated <- FALSE - - # Fleet - # Create the fishing fleet - fishing_fleet_selectivity <- new(fims$LogisticSelectivity) - fishing_fleet_selectivity$median$value <- om_input$sel_fleet$fleet1$A50.sel1 - fishing_fleet_selectivity$median$is_random_effect <- FALSE - fishing_fleet_selectivity$median$estimated <- TRUE - fishing_fleet_selectivity$slope$value <- om_input$sel_fleet$fleet1$slope.sel1 - fishing_fleet_selectivity$slope$is_random_effect <- FALSE - fishing_fleet_selectivity$slope$estimated <- TRUE - - fishing_fleet <- new(fims$Fleet) - fishing_fleet$nages <- om_input$nages - fishing_fleet$nyears <- om_input$nyr - fishing_fleet$log_Fmort <- log(om_output$f) - fishing_fleet$estimate_F <- TRUE - fishing_fleet$random_F <- FALSE - fishing_fleet$log_q <- log(1.0) - fishing_fleet$estimate_q <- FALSE - fishing_fleet$random_q <- FALSE - fishing_fleet$log_obs_error$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1))) - fishing_fleet$log_obs_error$estimated <- FALSE - # Need get_id() for setting up observed agecomp and index data? - fishing_fleet$SetAgeCompLikelihood(1) - fishing_fleet$SetIndexLikelihood(1) - fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id()) - fishing_fleet$SetObservedIndexData(fishing_fleet_index$get_id()) - fishing_fleet$SetObservedAgeCompData(fishing_fleet_age_comp$get_id()) - - # Create the survey fleet - survey_fleet_selectivity <- new(fims$LogisticSelectivity) - survey_fleet_selectivity$median$value <- om_input$sel_survey$survey1$A50.sel1 - survey_fleet_selectivity$median$is_random_effect <- FALSE - survey_fleet_selectivity$median$estimated <- TRUE - survey_fleet_selectivity$slope$value <- om_input$sel_survey$survey1$slope.sel1 - survey_fleet_selectivity$slope$is_random_effect <- FALSE - survey_fleet_selectivity$slope$estimated <- TRUE - - survey_fleet <- new(fims$Fleet) - survey_fleet$is_survey <- TRUE - survey_fleet$nages <- om_input$nages - survey_fleet$nyears <- om_input$nyr - # survey_fleet$log_Fmort <- rep(log(0.0000000000000000000000000001), om_input$nyr) #-Inf? - survey_fleet$estimate_F <- FALSE - survey_fleet$random_F <- FALSE - survey_fleet$log_q <- log(om_output$survey_q$survey1) - survey_fleet$estimate_q <- TRUE - survey_fleet$random_q <- FALSE - survey_fleet$log_obs_error$value <- log(sqrt(log(em_input$cv.survey$survey1^2 + 1))) - survey_fleet$log_obs_error$estimated <- FALSE - survey_fleet$SetAgeCompLikelihood(1) - survey_fleet$SetIndexLikelihood(1) - survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id()) - survey_fleet$SetObservedIndexData(survey_fleet_index$get_id()) - survey_fleet$SetObservedAgeCompData(survey_fleet_age_comp$get_id()) - - # Population - population <- new(fims$Population) - # is it a problem these are not Parameters in the Population interface? - # the Parameter class (from rcpp/rcpp_objects/rcpp_interface_base) cannot handle vectors, - # do we need a ParameterVector class? - population$log_M <- rep(log(om_input$M.age[1]), om_input$nyr * om_input$nages) - population$estimate_M <- FALSE - population$log_init_naa <- log(om_output$N.age[1, ]) - population$estimate_init_naa <- TRUE - population$nages <- om_input$nages - population$ages <- om_input$ages - population$nfleets <- sum(om_input$fleet_num, om_input$survey_num) - population$nseasons <- 1 - population$nyears <- om_input$nyr - population$prop_female <- om_input$proportion.female[1] - population$SetMaturity(maturity$get_id()) - population$SetGrowth(ewaa_growth$get_id()) - population$SetRecruitment(recruitment$get_id()) - - ## Set-up TMB - fims$CreateTMBModel() - parameters <- list(p = fims$get_fixed()) - obj<-TMB::MakeADFun(data = list(), parameters, DLL = "FIMS") - - fims$clear() - return(obj) +init_fims <- function(i) { + fims <- Rcpp::Module("fims", PACKAGE = "FIMS") + + # Recruitment + recruitment <- new(fims$BevertonHoltRecruitment) + # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is taken + # before the likelihood calculation + recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) + recruitment$log_rzero$value <- log(om_input$R0) # + runif(1,min=0, max=1000)) + recruitment$log_rzero$is_random_effect <- FALSE + recruitment$log_rzero$estimated <- TRUE + recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) + recruitment$logit_steep$is_random_effect <- FALSE + recruitment$logit_steep$estimated <- FALSE + recruitment$estimate_deviations <- TRUE + # recruit deviations should enter the model in normal space. + # The log is taken in the likelihood calculations + # alternative setting: recruitment$deviations <- rep(1, length(om_input$logR.resid)) + recruitment$deviations <- exp(om_input$logR.resid) + + # Data + catch <- em_input$L.obs$fleet1 + fishing_fleet_index <- new(fims$Index, length(catch)) + fishing_fleet_index$index_data <- catch + fishing_fleet_age_comp <- new(fims$AgeComp, length(catch), om_input$nages) + fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1 + + survey_index <- em_input$surveyB.obs$survey1 + survey_fleet_index <- new(fims$Index, length(survey_index)) + survey_fleet_index$index_data <- survey_index + survey_fleet_age_comp <- new(fims$AgeComp, length(survey_index), om_input$nages) + survey_fleet_age_comp$age_comp_data <- c(t(em_input$survey.age.obs$survey1)) * em_input$n.survey$survey1 + + # Growth + ewaa_growth <- new(fims$EWAAgrowth) + ewaa_growth$ages <- om_input$ages + ewaa_growth$weights <- om_input$W.mt + + # Maturity + maturity <- new(fims$LogisticMaturity) + maturity$median$value <- om_input$A50.mat + maturity$median$is_random_effect <- FALSE + maturity$median$estimated <- FALSE + maturity$slope$value <- om_input$slope + maturity$slope$is_random_effect <- FALSE + maturity$slope$estimated <- FALSE + + # Fleet + # Create the fishing fleet + fishing_fleet_selectivity <- new(fims$LogisticSelectivity) + fishing_fleet_selectivity$median$value <- om_input$sel_fleet$fleet1$A50.sel1 + fishing_fleet_selectivity$median$is_random_effect <- FALSE + fishing_fleet_selectivity$median$estimated <- TRUE + fishing_fleet_selectivity$slope$value <- om_input$sel_fleet$fleet1$slope.sel1 + fishing_fleet_selectivity$slope$is_random_effect <- FALSE + fishing_fleet_selectivity$slope$estimated <- TRUE + + fishing_fleet <- new(fims$Fleet) + fishing_fleet$nages <- om_input$nages + fishing_fleet$nyears <- om_input$nyr + fishing_fleet$log_Fmort <- log(om_output$f) + fishing_fleet$estimate_F <- TRUE + fishing_fleet$random_F <- FALSE + fishing_fleet$log_q <- log(1.0) + fishing_fleet$estimate_q <- FALSE + fishing_fleet$random_q <- FALSE + fishing_fleet$log_obs_error$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1))) + fishing_fleet$log_obs_error$estimated <- FALSE + # Need get_id() for setting up observed agecomp and index data? + fishing_fleet$SetAgeCompLikelihood(1) + fishing_fleet$SetIndexLikelihood(1) + fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id()) + fishing_fleet$SetObservedIndexData(fishing_fleet_index$get_id()) + fishing_fleet$SetObservedAgeCompData(fishing_fleet_age_comp$get_id()) + + # Create the survey fleet + survey_fleet_selectivity <- new(fims$LogisticSelectivity) + survey_fleet_selectivity$median$value <- om_input$sel_survey$survey1$A50.sel1 + survey_fleet_selectivity$median$is_random_effect <- FALSE + survey_fleet_selectivity$median$estimated <- TRUE + survey_fleet_selectivity$slope$value <- om_input$sel_survey$survey1$slope.sel1 + survey_fleet_selectivity$slope$is_random_effect <- FALSE + survey_fleet_selectivity$slope$estimated <- TRUE + + survey_fleet <- new(fims$Fleet) + survey_fleet$is_survey <- TRUE + survey_fleet$nages <- om_input$nages + survey_fleet$nyears <- om_input$nyr + # survey_fleet$log_Fmort <- rep(log(0.0000000000000000000000000001), om_input$nyr) #-Inf? + survey_fleet$estimate_F <- FALSE + survey_fleet$random_F <- FALSE + survey_fleet$log_q <- log(om_output$survey_q$survey1) + survey_fleet$estimate_q <- TRUE + survey_fleet$random_q <- FALSE + survey_fleet$log_obs_error$value <- log(sqrt(log(em_input$cv.survey$survey1^2 + 1))) + survey_fleet$log_obs_error$estimated <- FALSE + survey_fleet$SetAgeCompLikelihood(1) + survey_fleet$SetIndexLikelihood(1) + survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id()) + survey_fleet$SetObservedIndexData(survey_fleet_index$get_id()) + survey_fleet$SetObservedAgeCompData(survey_fleet_age_comp$get_id()) + + # Population + population <- new(fims$Population) + # is it a problem these are not Parameters in the Population interface? + # the Parameter class (from rcpp/rcpp_objects/rcpp_interface_base) cannot handle vectors, + # do we need a ParameterVector class? + population$log_M <- rep(log(om_input$M.age[1]), om_input$nyr * om_input$nages) + population$estimate_M <- FALSE + population$log_init_naa <- log(om_output$N.age[1, ]) + population$estimate_init_naa <- TRUE + population$nages <- om_input$nages + population$ages <- om_input$ages + population$nfleets <- sum(om_input$fleet_num, om_input$survey_num) + population$nseasons <- 1 + population$nyears <- om_input$nyr + population$prop_female <- om_input$proportion.female[1] + population$SetMaturity(maturity$get_id()) + population$SetGrowth(ewaa_growth$get_id()) + population$SetRecruitment(recruitment$get_id()) + + ## Set-up TMB + fims$CreateTMBModel() + parameters <- list(p = fims$get_fixed()) + obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS") + + fims$clear() + return(obj) } -run_fims<-function(id,begin,end){ - library(FIMS) - library(minimizR) - minimizer<-2 - results<-list() - - - index<-1 - for(i in begin[id]:end[id]){ - obj<-init_fims(i) - if(minimizer == 1){ - results[[index]] <- optim(obj$par, obj$fn, obj$gr, - method = "BFGS", - control = list(maxit = 1000000, reltol = 1e-15) - ) - }else{ - results[[index]]<-minimizR(obj$par, obj$fn, obj$gr, control = list(tolerance = 1e-8)) - } - - index = index+1 +run_fims <- function(id, begin, end) { + library(FIMS) + library(minimizR) + minimizer <- 2 + results <- list() + + + index <- 1 + for (i in begin[id]:end[id]) { + obj <- init_fims(i) + if (minimizer == 1) { + results[[index]] <- optim(obj$par, obj$fn, obj$gr, + method = "BFGS", + control = list(maxit = 1000000, reltol = 1e-15) + ) + } else { + results[[index]] <- minimizR(obj$par, obj$fn, obj$gr, control = list(tolerance = 1e-8)) } - return(results) + + index <- index + 1 + } + return(results) } @@ -203,32 +203,32 @@ nsims <- NUMBER_OF_MODEL_RUNS begin <- rep(0, ns) end <- rep(0, ns) -#create scenario segments +# create scenario segments if (id == 0) { - segments <- nsims / ns - print(paste("segments ", segments)) - for (i in 1:ns) { - if (i < ns) { - begin[i] <- as.integer((i - 1) * segments+1) - end[i] <- as.integer(i * segments) - } else{ - begin[i] <- as.integer((i - 1) * segments+1) - end[i] <- nsims - } + segments <- nsims / ns + print(paste("segments ", segments)) + for (i in 1:ns) { + if (i < ns) { + begin[i] <- as.integer((i - 1) * segments + 1) + end[i] <- as.integer(i * segments) + } else { + begin[i] <- as.integer((i - 1) * segments + 1) + end[i] <- nsims } - print(begin) - print(end) + } + print(begin) + print(end) } -completed<-list() +completed <- list() -sfInit(parallel=TRUE, cpus=ns) -sfExport( "run_fims", "init_fims", "om_input","om_output", "em_input") -start<-Sys.time() +sfInit(parallel = TRUE, cpus = ns) +sfExport("run_fims", "init_fims", "om_input", "om_output", "em_input") +start <- Sys.time() append(sfLapply(1:ns, run_fims, begin, end), completed) -end<-Sys.time() -runtime<-end-start -print(paste0(paste0(paste0(NUMBER_OF_MODEL_RUNS," model runs completed in "),runtime)," seconds.")) +end <- Sys.time() +runtime <- end - start +print(paste0(paste0(paste0(NUMBER_OF_MODEL_RUNS, " model runs completed in "), runtime), " seconds.")) -line=paste("Snowfall runtime ", runtime) -write(line,file="time.txt",append=TRUE) +line <- paste("Snowfall runtime ", runtime) +write(line, file = "time.txt", append = TRUE) diff --git a/tests/testthat/test-data-object.R b/tests/testthat/test-data-object.R index 66532750..d3be61c5 100644 --- a/tests/testthat/test-data-object.R +++ b/tests/testthat/test-data-object.R @@ -23,7 +23,6 @@ fleet_names_index <- dplyr::filter( nindex <- length(fleet_names_index) test_that("Can add index data to model", { - indexdat <- vector(mode = "list", length = nindex) names(indexdat) <- fleet_names_index @@ -38,7 +37,6 @@ test_that("Can add index data to model", { }) test_that("Can add agecomp data to model", { - agecompdat <- vector(mode = "list", length = nagecomp) names(agecompdat) <- fleet_names_agecomp diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index d79a02dc..d7bfaa1c 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -16,7 +16,7 @@ FIMS_C0_estimation <- ASSAMC::save_initial_input( case_name = "FIMS_C0_estimation" ) -# generate om_input, om_output, and em_input +# generate om_input, om_output, and em_input # using function from the model comparison project ASSAMC::run_om(input_list = FIMS_C0_estimation) @@ -37,10 +37,10 @@ setup_fims <- function(om_input, om_output, em_input) { # when there are other options, this would be where the option would be chosen) test_env$recruitment <- methods::new(test_env$fims$BevertonHoltRecruitment) - # NOTE: in first set of parameters below (for recruitment), - # $is_random_effect (default is FALSE) and $estimated (default is FALSE) - # are defined even if they match the defaults in order to provide an example - # of how that is done. Other sections of the code below leave defaults in + # NOTE: in first set of parameters below (for recruitment), + # $is_random_effect (default is FALSE) and $estimated (default is FALSE) + # are defined even if they match the defaults in order to provide an example + # of how that is done. Other sections of the code below leave defaults in # place as appropriate. # set up logR_sd @@ -186,14 +186,14 @@ test_that("deterministic test of fims", { om_output = om_output, em_input = em_input ) - # Set-up + # Set-up deterministic_env$fims$CreateTMBModel() # CreateTMBModel calls a function in information that loops # over all the populations and fleets and sets all the pointers parameters <- list(p = deterministic_env$fims$get_fixed()) # get_fixed function is an Rcpp function that loops over all Rcpp # modules and returned a vector of parameters being estimated - + # Set up TMB's computational graph obj <- MakeADFun(data = list(), parameters, DLL = "FIMS") @@ -354,7 +354,7 @@ test_that("nll test of fims", { # recruitment likelihood # log_devs is of length nyr-1 rec_nll <- -sum(dnorm( - nll_env$recruitment$log_devs, rep(0, om_input$nyr-1), + nll_env$recruitment$log_devs, rep(0, om_input$nyr - 1), om_input$logR_sd, TRUE )) @@ -485,15 +485,15 @@ test_that("estimation test of fims", { # recruitment log deviations # the initial value of om_input$logR.resid is dropped from the model sdr_rdev <- sdr_report[which(rownames(sdr_report) == "LogRecDev"), ] - rdev_are <- rep(0, length(om_input$logR.resid)-1) + rdev_are <- rep(0, length(om_input$logR.resid) - 1) - for (i in 1:length(report$log_recruit_dev[[1]]) ){ - rdev_are[i] <- abs(report$log_recruit_dev[[1]][i] - om_input$logR.resid[i+1]) # / + for (i in 1:length(report$log_recruit_dev[[1]])) { + rdev_are[i] <- abs(report$log_recruit_dev[[1]][i] - om_input$logR.resid[i + 1]) # / # exp(om_input$logR.resid[i]) # expect_lte(rdev_are[i], 1) # 1 } expect_lte( - sum(rdev_are > qnorm(.975) * sdr_rdev[1:length(om_input$logR.resid)-1, 2]), + sum(rdev_are > qnorm(.975) * sdr_rdev[1:length(om_input$logR.resid) - 1, 2]), 0.05 * length(om_input$logR.resid) ) @@ -622,14 +622,14 @@ test_that("run FIMS in a for loop", { # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is taken # before the likelihood calculation recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) - recruitment$log_rzero$value <- 13 #log(om_input$R0) + recruitment$log_rzero$value <- 13 # log(om_input$R0) recruitment$log_rzero$is_random_effect <- FALSE recruitment$log_rzero$estimated <- TRUE recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) recruitment$logit_steep$is_random_effect <- FALSE recruitment$logit_steep$estimated <- FALSE recruitment$estimate_log_devs <- TRUE - recruitment$log_devs <- rep(0, length(om_input$logR.resid)-1) + recruitment$log_devs <- rep(0, length(om_input$logR.resid) - 1) # Data catch <- em_input$L.obs$fleet1 @@ -736,12 +736,12 @@ test_that("run FIMS in a for loop", { fims$CreateTMBModel() parameters <- list(p = fims$get_fixed()) obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS") - + opt <- with(obj, optim(par, fn, gr, method = "BFGS", control = list(maxit = 1000000, reltol = 1e-15) )) - + report <- obj$report(obj$par) expect_false(is.null(report)) diff --git a/tests/testthat/test-fims-rcpp.R b/tests/testthat/test-fims-rcpp.R index ca219bea..58cd924a 100644 --- a/tests/testthat/test-fims-rcpp.R +++ b/tests/testthat/test-fims-rcpp.R @@ -1,5 +1,4 @@ test_that("Rcpp interface works for modules", { - expect_no_error(parameter <- new(Parameter, .1)) expect_no_error(beverton_holt <- new(BevertonHoltRecruitment)) expect_no_error(logistic_selectivity <- new(LogisticSelectivity)) diff --git a/tests/testthat/test-fleet-interface.R b/tests/testthat/test-fleet-interface.R index 0562e5db..637b29a8 100644 --- a/tests/testthat/test-fleet-interface.R +++ b/tests/testthat/test-fleet-interface.R @@ -1,6 +1,5 @@ test_that("Fleet: selectivity IDs can be added to the fleet module", { - # Create selectivity for fleet 1 selectivity_fleet1 <- new(LogisticSelectivity) expect_equal((selectivity_fleet1$get_id()), 1) diff --git a/tests/testthat/test-maturity-interface.R b/tests/testthat/test-maturity-interface.R index 907ef968..4ded4336 100644 --- a/tests/testthat/test-maturity-interface.R +++ b/tests/testthat/test-maturity-interface.R @@ -1,5 +1,4 @@ test_that("Maturity input settings work as expected", { - # Create maturity1 maturity1 <- new(LogisticMaturity) diff --git a/tests/testthat/test-recruitment-interface.R b/tests/testthat/test-recruitment-interface.R index 4d55d639..55e122b8 100644 --- a/tests/testthat/test-recruitment-interface.R +++ b/tests/testthat/test-recruitment-interface.R @@ -1,6 +1,5 @@ library(testthat) test_that("Recruitment input settings work as expected", { - # Create recruitment recruitment <- new(BevertonHoltRecruitment) h <- 0.75 diff --git a/tests/testthat/test-selectivity-interface.R b/tests/testthat/test-selectivity-interface.R index d720a0ed..398fe0c9 100644 --- a/tests/testthat/test-selectivity-interface.R +++ b/tests/testthat/test-selectivity-interface.R @@ -1,5 +1,4 @@ test_that("Selectivity input settings work as expected", { - # Create selectivity1 selectivity1 <- new(LogisticSelectivity) diff --git a/vignettes/fims-demo.Rmd b/vignettes/fims-demo.Rmd index 5053674f..8aab38bc 100644 --- a/vignettes/fims-demo.Rmd +++ b/vignettes/fims-demo.Rmd @@ -257,13 +257,13 @@ We also need to set up log recruitment deviations. FIMS recruitment modules have ```{r rec-devs} recruitment$estimate_log_devs <- FALSE recruitment$log_devs <- c( - 0.08904850, 0.43787763, -0.13299042, -0.43251973, - 0.64861200, 0.50640852, -0.06958319, 0.30246260, - -0.08257384, 0.20740372, 0.15289604, -0.21709207, - -0.13320626, 0.11225374, -0.10650836, 0.26877132, - 0.24094126, -0.54480751, -0.23680557, -0.58483386, - 0.30122785, 0.21930545, -0.22281699, -0.51358369, - 0.15740234, -0.53988240, -0.19556523, 0.20094360, + 0.08904850, 0.43787763, -0.13299042, -0.43251973, + 0.64861200, 0.50640852, -0.06958319, 0.30246260, + -0.08257384, 0.20740372, 0.15289604, -0.21709207, + -0.13320626, 0.11225374, -0.10650836, 0.26877132, + 0.24094126, -0.54480751, -0.23680557, -0.58483386, + 0.30122785, 0.21930545, -0.22281699, -0.51358369, + 0.15740234, -0.53988240, -0.19556523, 0.20094360, 0.37248740, -0.07163145 ) ```