Add list form for pollock and some tests -- matches but crashes often

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)
+
+
+

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)