diff --git a/content/R/pk_prepare_FIMS_inputs_by_process.R b/content/R/pk_prepare_FIMS_inputs_by_process.R new file mode 100644 index 0000000..d268b6e --- /dev/null +++ b/content/R/pk_prepare_FIMS_inputs_by_process.R @@ -0,0 +1,179 @@ + + +## define the dimensions and global variables + + +setup_fleets <- function(fleets,population, catchability, selectivity, data){ + flts <- list() + for(i in 1:length(fleets$names)){ + isfleet <- fleets$type[i] == 'landings' + flts[[i]] <- methods::new(Fleet) + flts[[i]]$nages <- length(population$ages) + flts[[i]]$nyears <- population$nyears + fltindex <- methods::new(Index, nyears) + fltacomp <- methods::new(AgeComp, flts[[i]]$nyears, flts[[i]]$nages) + fltdat <- filter(data@data, name==fleets$names[i]) + acomp <- FIMS::m_agecomp(age_frame, fleets$names[i]) + if(isfleet){ + flts[[i]]$is_survey <- FALSE + flts[[i]]$log_Fmort <- log(fleets$init_Fmort[[i]]) + flts[[i]]$estimate_F <- fleets$estimate_F[i] + flts[[i]]$random_F <- fleets$reFmort[i] + fltindex$index_data <- FIMS::m_landings(data) + flts[[i]]$log_obs_error <- log(filter(fltdat, type=='landings')$uncertainty) + fltacomp$age_comp_data <- acomp * filter(fltdat, type=='age')$uncertainty + } else { + flts[[i]]$is_survey <- TRUE + flts[[i]]$log_q <- catchability$init_pars[i] + flts[[i]]$estimate_q <- is.na(catchability$fix_pars[i]) + flts[[i]]$random_q <- catchability$re[i]!='none' + fltindex$index_data <- FIMS::m_index(data, fleets$names[i]) + flts[[i]]$log_obs_error <- log(filter(fltdat, type=='index')$uncertainty) + fltacomp$age_comp_data <- acomp * filter(fltdat, type=='age')$uncertainty + } + ## Set Index, AgeComp, and Selectivity using the IDs from the modules defined above + flts[[i]]$SetObservedIndexData(fltindex$get_id()) + flts[[i]]$SetObservedAgeCompData(fltacomp$get_id()) + flts[[i]]$SetSelectivity(selectivity[[i]]$get_id()) + } + message("Finished setting up fleets") +} + +setup_selex <- function(selectivity){ + s <- selectivity + sels <- list() + for(i in 1:length(s$names)){ + if(s$form[i]=='DoubleLogistic'){ + sels[[i]] <- methods::new(DoubleLogisticSelectivity) + + sels[[i]]$inflection_point_asc$value <- s$init_pars[[i]][1] + sels[[i]]$inflection_point_asc$is_random_effect <- s$re[1]!='none' & 1 %in% s$where[[i]] + sels[[i]]$inflection_point_asc$estimated <- !(1 %in% s$fix_pars[[i]]) + + sels[[i]]$slope_asc$value <- s$init_pars[[i]][2] + sels[[i]]$slope_asc$is_random_effect <- s$re[1]!='none' & 2 %in% s$where[[i]] + sels[[i]]$slope_asc$estimated <- !(2 %in% s$fix_pars[[i]]) + + sels[[i]]$inflection_point_desc$value <- s$init_pars[[i]][3] + sels[[i]]$inflection_point_desc$is_random_effect <- s$re[1]!='none' & 3 %in% s$where[[i]] + sels[[i]]$inflection_point_desc$estimated <- !(1 %in% s$fix_pars[[i]]) + + sels[[i]]$slope_desc$value <- s$init_pars[[i]][4] + sels[[i]]$slope_desc$is_random_effect <- s$re[1]!='none' & 4 %in% s$where[[i]] + sels[[i]]$slope_desc$estimated <- !(4 %in% s$fix_pars[[i]]) + } else if(s$form[i]=='Logistic'){ + sels[[i]] <- methods::new(LogisticSelectivity) + + sels[[i]]$inflection_point$value <- s$init_pars[[i]][1] + sels[[i]]$inflection_point$is_random_effect <- s$re[1]!='none' & 1 %in% s$where[[i]] + sels[[i]]$inflection_point$estimated <- !(1 %in% s$fix_pars[[i]]) + + sels[[i]]$slope$value <- s$init_pars[[i]][2] + sels[[i]]$slope$is_random_effect <- s$re[1]!='none' & 2 %in% s$where[[i]] + sels[[i]]$slope$estimated <- !(2 %in% s$fix_pars[[i]]) + } + } + message("Finished setting up selex") + return(sels) +} +setup_recruitment <- function(recruitment){ + rec <- recruitment + rec <- methods::new(BevertonHoltRecruitment) + ## methods::show(BevertonHoltRecruitment) + rec$log_sigma_recruit$value <- log(recruitment$init_pars[3]) + rec$log_rzero$value <- log(recruitment$init_pars[1]) + rec$log_rzero$is_random_effect <- FALSE + rec$log_rzero$estimated <- (1 %in% recruitment$fix_pars) + steep <- recruitment$init_pars[2] + rec$logit_steep$value <- -log(1.0 - steep) + log(steep - 0.2) + rec$logit_steep$is_random_effect <- FALSE + rec$logit_steep$estimated <- (2 %in% recruitment$fix_pars) + rec$estimate_log_devs <- estimate_recdevs + rec$log_devs <- recruitment$init_recdevs + message("Finished setting up recruitment") + return(rec) +} +setup_growth <- function(growth){ + gr <- methods::new(EWAAgrowth) + gr$ages <- ages + gr$weights <- growth$waa + message("Finished setting up growth") + return(gr) +} +setup_maturity <- function(maturity){ + mat <- new(LogisticMaturity) + mat$inflection_point$value <- maturity$init_pars[1] + mat$inflection_point$is_random_effect <- FALSE + mat$inflection_point$estimated <- !(1 %in% maturity$fix_pars) + mat$slope$value <- maturity$init_pars[2] + mat$slope$is_random_effect <- FALSE + mat$slope$estimated <- !(2 %in% maturity$fix_pars) + message("Finished setting up maturity") + return(mat) +} +setup_population <- function(population, matout, growthout, recout){ + pop <- new(Population) + pop$log_M <- log(population$M) + pop$estimate_M <- FALSE + pop$log_init_naa <- log(population$init_naa) + pop$estimate_init_naa <- population$fix_naa + pop$nages <- population$nages + pop$ages <- population$ages + pop$nfleets <- population$nfleets + pop$nseasons <- 1 + pop$nyears <- population$nyears + ## pop$prop_female <- 1.0 # ASAP assumption + pop$SetMaturity(matout$get_id()) + pop$SetGrowth(growthout$get_id()) + pop$SetRecruitment(recout$get_id()) + message("Finished setting up population and linking maturity, growth, recruitment") +} + +## Order of parameters for DL is asc inf, asc slop, desc inf, +## desc slope +mle <- parfinal +data <- age_frame +myNA <- rep(NA,3) +population <- list(ages=ages, nfleets=4, nyears=nyears, + years=years, nages=length(ages), + M=as.numeric(t(matrix(rep(pkfitfinal$rep$M, each=nyears), nrow=nyears))), + init_naa= exp(c(log(pkfitfinal$rep$recruit[1]), log(pkfitfinal$rep$initN)) + log(1e9)), + fix_naa=TRUE) +fleets <- list(names=c('fleet1', 'survey2', 'survey3', 'survey6'), + type=c('landings', rep('index',3)), + estimate_Fmort=c(TRUE, myNA), + init_Fmort=list(pkfitfinal$rep$F, myNA), + reFmort=c(FALSE, myNA)) +sel <- list(names=fleets$names, + form=c('DoubleLogistic', 'DoubleLogistic', 'Logistic', 'DoubleLogistic'), + fix_pars=list(NA, 3:4, NA, 1:2), + re=rep('none',4), + init_pars=list(c(mle$inf1_fsh_mean, exp(mle$log_slp1_fsh_mean), mle$inf2_fsh_mean, exp(mle$log_slp2_fsh_mean)), + c(mle$inf1_srv2, exp(mle$log_slp1_srv2), mle$inf2_srv2, exp(mle$log_slp2_srv2)), + c(mle$inf1_srv3, exp(mle$log_slp1_srv3)), + c(mle$inf1_srv6, exp(mle$log_slp1_srv6), mle$inf2_srv6, exp(mle$log_slp2_srv6)))) +catchability <- list(names=fleets$names, + form=c('none', rep('linear', times=3)), + re=rep('none', 4), + init_pars=c(NA, mle$log_q2_mean, mle$log_q3_mean, mle$log_q6), + ## init_repars=c('rho'=.8, 'sd'=.1), + fix_pars=c(NA,myNA), + fix_repars=c(NA,myNA)) +recruitment <- list(type='BevertonHolt', + ## R0, steepness, sigmaR + init_pars=c(exp(parfinal$mean_log_recruit + log(1e9)), .99999, parfinal$sigmaR), + fix_pars=c(2:3), + init_recdevs=parfinal$dev_log_recruit[-1], + re=c('none', 'none')) +growth <- list(type='EWAA', waa=pkinput$dat$wt_srv1[1,]) +maturity <- list(type='logistic', init_pars=c(4.5, 1.5), fix_pars=1:2) + +selexout <- setup_selex(sel) +setup_fleets(fleets, population, catchability, selectivity=selexout, data=age_frame) +recout <- setup_recruitment(recruitment) +growthout <- setup_growth(growth) +matout <- setup_maturity(maturity) +setup_population(population, matout, growthout, recout) + + + diff --git a/content/R/pk_prepare_dat.R b/content/R/pk_prepare_dat.R new file mode 100644 index 0000000..e680b99 --- /dev/null +++ b/content/R/pk_prepare_dat.R @@ -0,0 +1,179 @@ + + +## build a FIMS and PK data set that match + +pkfitfinal <- readRDS("data_files/pkfitfinal.RDS") +pkfit0 <- readRDS("data_files/pkfit0.RDS") +parfinal <- pkfitfinal$obj$env$parList() +pkinput0 <- readRDS('data_files/pkinput0.RDS') +fimsdat <- pkdat0 <- pkinput0$dat +pkinput <- readRDS('data_files/pkinput.RDS') + + + +## need to fill missing years with -999 so it's ignored in FIMS +ind2 <- 0 * pkfit0$rep$Eindxsurv2 - 999 +ind2[which(years %in% fimsdat$srvyrs2)] <- fimsdat$indxsurv2 +CV2 <- rep(1, length = nyears) # actually SE in log space +CV2[which(years %in% fimsdat$srvyrs2)] <- fimsdat$indxsurv_log_sd2 +paa2 <- pkfit0$rep$Esrvp2 * 0 - 999 +paa2[which(years %in% fimsdat$srv_acyrs2), ] <- fimsdat$srvp2 +Npaa2 <- rep(1, nyears) +Npaa2[which(years %in% fimsdat$srv_acyrs2)] <- fimsdat$multN_srv2 + +ind3 <- 0 * pkfit0$rep$Eindxsurv3 - 999 +ind3[which(years %in% fimsdat$srvyrs3)] <- fimsdat$indxsurv3 +CV3 <- rep(1, length = nyears) # actually SE in log space +CV3[which(years %in% fimsdat$srvyrs3)] <- fimsdat$indxsurv_log_sd3 +paa3 <- pkfit0$rep$Esrvp3 * 0 - 999 +paa3[which(years %in% fimsdat$srv_acyrs3), ] <- fimsdat$srvp3 +Npaa3 <- rep(1, nyears) +Npaa3[which(years %in% fimsdat$srv_acyrs3)] <- fimsdat$multN_srv3 + +ind6 <- 0 * pkfit0$rep$Eindxsurv6 - 999 +ind6[which(years %in% fimsdat$srvyrs6)] <- fimsdat$indxsurv6 +CV6 <- rep(1, length = nyears) # actually SE in log space +CV6[which(years %in% fimsdat$srvyrs6)] <- fimsdat$indxsurv_log_sd6 +paa6 <- pkfit0$rep$Esrvp6 * 0 - 999 +paa6[which(years %in% fimsdat$srv_acyrs6), ] <- fimsdat$srvp6 +Npaa6 <- rep(1, nyears) +Npaa6[which(years %in% fimsdat$srv_acyrs6)] <- fimsdat$multN_srv6 + +## repeat with fish catch at age, using expected in missing years +caa <- pkfit0$rep$Ecatp * 0 - 999 +caa[which(years %in% fimsdat$fshyrs), ] <- fimsdat$catp +Ncaa <- rep(1, nyears) +Ncaa[which(years %in% fimsdat$fshyrs)] <- fimsdat$multN_fsh + + + +## put into fims friendly form +res <- data.frame( + type = character(), + name = character(), + age = integer(), + datestart = character(), + dateend = character(), + value = double(), + unit = character(), + uncertainty = double() +) +landings <- data.frame( + type = "landings", + name = "fleet1", + age = NA, + datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"), + dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"), + value = as.numeric(fimsdat$cattot) * 1e3, + unit = "mt", + uncertainty = fimsdat$cattot_log_sd[1] +) +index2 <- data.frame( + type = "index", + name = "survey2", + age = NA, + datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"), + dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"), + value = ifelse(ind2 > 0, ind2 * 1e9, ind2), + unit = "", + uncertainty = CV2 +) +index3 <- data.frame( + type = "index", + name = "survey3", + age = NA, + datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"), + dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"), + value = ifelse(ind3 > 0, ind3 * 1e9, ind3), + unit = "", + uncertainty = CV3 +) +index6 <- data.frame( + type = "index", + name = "survey6", + age = NA, + datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"), + dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"), + value = ifelse(ind6 > 0, ind6 * 1e9, ind6), + unit = "", + uncertainty = CV6 +) +## these have -999 for missing data years +catchage <- data.frame( + type = "age", + name = "fleet1", + age = rep(seq(1, nages), nyears), + datestart = rep(paste0( + seq(fimsdat$styr, fimsdat$endyr), "-01-01" + ), each = nages), + dateend = rep(paste0( + seq(fimsdat$styr, fimsdat$endyr), "-12-31" + ), each = nages), + value = as.numeric(t(caa)), + unit = "", + uncertainty = rep(Ncaa, each = nages) +) +indexage2 <- data.frame( + type = "age", + name = "survey2", + age = rep(seq(1, nages), nyears), + datestart = rep(paste0( + seq(fimsdat$styr, fimsdat$endyr), "-01-01" + ), each = nages), + dateend = rep(paste0( + seq(fimsdat$styr, fimsdat$endyr), "-12-31" + ), each = nages), + value = as.numeric(t(paa2)), + unit = "", + uncertainty = rep(Npaa2, each = nages) +) +indexage3 <- data.frame( + type = "age", + name = "survey3", + age = rep(seq(1, nages), nyears), + datestart = rep(paste0( + seq(fimsdat$styr, fimsdat$endyr), "-01-01" + ), each = nages), + dateend = rep(paste0( + seq(fimsdat$styr, fimsdat$endyr), "-12-31" + ), each = nages), + value = as.numeric(t(paa3)), + unit = "", + uncertainty = rep(Npaa3, each = nages) +) +indexage6 <- data.frame( + type = "age", + name = "survey6", + age = rep(seq(1, nages), nyears), + datestart = rep(paste0( + seq(fimsdat$styr, fimsdat$endyr), "-01-01" + ), each = nages), + dateend = rep(paste0( + seq(fimsdat$styr, fimsdat$endyr), "-12-31" + ), each = nages), + value = as.numeric(t(paa6)), + unit = "", + uncertainty = rep(Npaa6, each = nages) +) +indexage <- rbind(indexage2, indexage3, indexage6) +index <- rbind(index2, index3, index6) +## indexage=indexage2 +## index=index2 +res <- rbind(res, landings, index, catchage, indexage) +## rm(landings, index, catchage, indexage) + +age_frame <- FIMS::FIMSFrame(res) +fishery_catch <- FIMS::m_landings(age_frame) +fishery_agecomp <- FIMS::m_agecomp(age_frame, "fleet1") +survey_index2 <- FIMS::m_index(age_frame, "survey2") +survey_agecomp2 <- FIMS::m_agecomp(age_frame, "survey2") +survey_index3 <- FIMS::m_index(age_frame, "survey3") +survey_agecomp3 <- FIMS::m_agecomp(age_frame, "survey3") +survey_index6 <- FIMS::m_index(age_frame, "survey6") +survey_agecomp6 <- FIMS::m_agecomp(age_frame, "survey6") +# need to think about how to deal with multiple fleets - only using 1 fleeet for now +fish_index <- methods::new(Index, nyears) +fish_age_comp <- methods::new(AgeComp, nyears, nages) +fish_index$index_data <- fishery_catch +fish_age_comp$age_comp_data <- + fishery_agecomp * catchage$uncertainty#rep(Ncaa, each=nages) diff --git a/content/run_pollock_tests.R b/content/run_pollock_tests.R new file mode 100644 index 0000000..4871ef0 --- /dev/null +++ b/content/run_pollock_tests.R @@ -0,0 +1,50 @@ + + +packages <- c("dplyr", "tidyr", "ggplot2", "TMB", "reshape2", "here", "remotes", "lubridate") +invisible(lapply(packages, library, character.only = TRUE)) +## install.packages("FIMS", repos = c("https://noaa-fims.r-universe.dev", "https://cloud.r-project.org")) +library(FIMS) +years <- 1970:2023 +nyears <- length(years) +nseasons <- 1 +nages <- 10 +ages <- 1:nages +## nfleets <- 1 +## This will fit the models bridging to FIMS (simplifying) +## source("fit_bridge_models.R") +## compare changes to model +source("R/pk_prepare_dat.R") +## Some global settings which I Think we can ignore for now +estimate_fish_selex <- TRUE +estimate_survey_selex <- TRUE +estimate_q2 <- TRUE +estimate_q3 <- TRUE +estimate_q6 <- TRUE +estimate_F <- TRUE +estimate_recdevs <- TRUE + + + +source("R/pk_prepare_FIMS_inputs.R") +## make FIMS model +success <- CreateTMBModel() +parameters <- list(p = get_fixed()) +obj1 <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE) +rep1 <- obj1$report() + +## Organize input lists by process +clear() +clear_logs() +source("R/pk_prepare_FIMS_inputs_by_process.R") +## make FIMS model +success <- CreateTMBModel() +parameters <- list(p = get_fixed()) +obj2 <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE) +## report values for the two models +rep2 <- obj2$report() + +## Test that the two models are identical +all.equal(rep1,rep2) + +clear() +clear_logs()