From 12abdefb22e7b8b5f0640b0389b5e07fc7a00877 Mon Sep 17 00:00:00 2001 From: brianstock-NOAA Date: Wed, 25 Aug 2021 17:56:03 -0400 Subject: [PATCH] fix errors in prepare_wham_input when missing indices or paa (butterfish, winter flounder) --- R/prepare_wham_input.R | 133 ++++++++++++++++++++--------------------- R/set_indices.R | 55 +++++++++-------- R/set_selectivity.R | 12 ++-- 3 files changed, 101 insertions(+), 99 deletions(-) diff --git a/R/prepare_wham_input.R b/R/prepare_wham_input.R index f518c983..1ebbcc5d 100644 --- a/R/prepare_wham_input.R +++ b/R/prepare_wham_input.R @@ -275,51 +275,51 @@ prepare_wham_input <- function(asap3 = NULL, model_name="WHAM for unnamed stock" par = list() map = list() random = character() - input = list( - data = data, - par = par, - map = map, - random = random, - years = NULL, years_full = NULL, ages.lab = NULL, model_name = model_name, asap3 = asap3) + input = list( + data = data, + par = par, + map = map, + random = random, + years = NULL, years_full = NULL, ages.lab = NULL, model_name = model_name, asap3 = asap3) - if(is.null(basic_info)) basic_info = list(recruit_model = recruit_model) - else basic_info$recruit_model = recruit_model + if(is.null(basic_info)) basic_info = list(recruit_model = recruit_model) + else basic_info$recruit_model = recruit_model - waa_opts = NULL - waa_names = c("waa") - if(any(names(basic_info) %in% waa_names)) waa_opts = basic_info[waa_names] + waa_opts = NULL + waa_names = c("waa") + if(any(names(basic_info) %in% waa_names)) waa_opts = basic_info[waa_names] catch_opts = NULL - catch_names = c("n_fleets","agg_catch", "catch_paa", "catch_cv","catch_Neff", "use_catch_paa", "selblock_pointer_fleets") - if(any(names(basic_info) %in% catch_names)) catch_opts = basic_info[catch_names] + catch_names = c("n_fleets","agg_catch", "catch_paa", "catch_cv","catch_Neff", "use_catch_paa", "selblock_pointer_fleets") + if(any(names(basic_info) %in% catch_names)) catch_opts = basic_info[catch_names] - index_opts = NULL - index_names = c("n_indices", "agg_indices", "agg_index_paa", "fracyr_indices", "index_cv", "index_Neff", "units_indices", - "units_index_paa", "use_indices", "use_index_paa", "selblock_pointer_indices") - if(any(names(basic_info) %in% index_names)) index_opts = basic_info[index_names] + index_opts = NULL + index_names = c("n_indices", "agg_indices", "agg_index_paa", "fracyr_indices", "index_cv", "index_Neff", "units_indices", + "units_index_paa", "use_indices", "use_index_paa", "selblock_pointer_indices") + if(any(names(basic_info) %in% index_names)) index_opts = basic_info[index_names] - F_opts = NULL - F_names = c("F") - if(any(names(basic_info) %in% F_names)) F_opts = basic_info[F_names] + F_opts = NULL + F_names = c("F") + if(any(names(basic_info) %in% F_names)) F_opts = basic_info[F_names] - waa_opts = NULL - waa_names = ("waa") - if(any(names(basic_info) %in% waa_names)) waa_opts = basic_info[waa_names] + waa_opts = NULL + waa_names = ("waa") + if(any(names(basic_info) %in% waa_names)) waa_opts = basic_info[waa_names] - q_opts = NULL - q_names = c("q","q_lower","q_upper") - if(any(names(basic_info) %in% q_names)) q_opts = basic_info[q_names] + q_opts = NULL + q_names = c("q","q_lower","q_upper") + if(any(names(basic_info) %in% q_names)) q_opts = basic_info[q_names] if(!is.null(asap3)) { asap3 = asap3$dat - input$asap3 = asap3 + input$asap3 = asap3 input$data$n_ages = asap3$n_ages input$data$fracyr_SSB = rep(asap3$fracyr_spawn, asap3$n_years) input$data$mature = asap3$maturity input$data$Fbar_ages = seq(asap3$Frep_ages[1], asap3$Frep_ages[2]) - input$years <- asap3$year1 + 1:asap3$n_years - 1 + input$years <- asap3$year1 + 1:asap3$n_years - 1 } else { @@ -327,63 +327,62 @@ prepare_wham_input <- function(asap3 = NULL, model_name="WHAM for unnamed stock" input = add_basic_info(input, basic_info) } - # print("start") + # print("start") #some basic input elements see the function code below input = initial_input_fn(input, basic_info) - # Catch - input = set_catch(input, catch_opts) - #print("catch") + # Catch + input = set_catch(input, catch_opts) + #print("catch") - # Indices/surveys - input = set_indices(input, index_opts) - #print("indices") + # Indices/surveys + input = set_indices(input, index_opts) + #print("indices") - # WAA in case we want to modify how weight-at age is handled - input = set_WAA(input, waa_opts) - #print("WAA") + # WAA in case we want to modify how weight-at age is handled + input = set_WAA(input, waa_opts) + #print("WAA") - # NAA and recruitment options - input = set_NAA(input, NAA_re) - #print("NAA") + # NAA and recruitment options + input = set_NAA(input, NAA_re) + #print("NAA") - input = set_q(input, q_opts) - #print("q") + input = set_q(input, q_opts) + #print("q") - # Selectivity - input = set_selectivity(input, selectivity) - #print("selectivity") + # Selectivity + input = set_selectivity(input, selectivity) + #print("selectivity") - # Age composition model - input = set_age_comp(input, age_comp) - #print("age_comp") + # Age composition model + input = set_age_comp(input, age_comp) + #print("age_comp") - #in case we want to add alternative F options - input = set_F(input, F_opts) - #print("F") + #in case we want to add alternative F options + input = set_F(input, F_opts) + #print("F") - #set up natural mortality - input = set_M(input, M) - #print("M") + #set up natural mortality + input = set_M(input, M) + #print("M") #set up ecov data and parameters. Probably want to make sure to do this after set_NAA. input = set_ecov(input, ecov) - #print("ecov") + #print("ecov") - # add vector of all observations for one step ahead residuals ========================== - input = set_osa_obs(input) - #print("osa_obs") + # add vector of all observations for one step ahead residuals ========================== + input = set_osa_obs(input) + #print("osa_obs") + # projection data will always be modified by 'prepare_projection' + input = set_proj(input, proj.opts = NULL) #proj options are used later after model fit, right? + #print("proj") - # projection data will always be modified by 'prepare_projection' - input = set_proj(input, proj.opts = NULL) #proj options are used later after model fit, right? - #print("proj") + #set any parameters as random effects + input = set_map(input) + #print("map") - #set any parameters as random effects - input = set_map(input) - #print("map") - - return(input) + return(input) #return(list(data=data, par = par, map = map, random = random, years = model_years, years_full = model_years, # ages.lab = paste0(1:data$n_ages, c(rep("",data$n_ages-1),"+")), model_name = model_name)) } diff --git a/R/set_indices.R b/R/set_indices.R index 2230cacb..259be5e4 100644 --- a/R/set_indices.R +++ b/R/set_indices.R @@ -3,36 +3,36 @@ set_indices = function(input, index_opts=NULL) data = input$data if(is.null(input$asap3)) { asap3 = NULL - if(is.null(index_opts$n_indices)) data$n_indices = 1 - else data$n_indices = index_opts$n_indices + if(is.null(index_opts$n_indices)) data$n_indices = 1 + else data$n_indices = index_opts$n_indices } - else { - asap3 = input$asap3 - which_indices <- which(asap3$use_index ==1) - asap3$n_indices = length(which_indices) - asap3$survey_index_units <- asap3$index_units[which_indices] - asap3$survey_acomp_units <- asap3$index_acomp_units[which_indices] - asap3$survey_WAA_pointers <- asap3$index_WAA_pointers[which_indices] - asap3$survey_month <- matrix(asap3$index_month[which_indices], asap3$n_years, asap3$n_indices, byrow = TRUE) - asap3$use_survey_acomp <- asap3$use_index_acomp[which_indices] - asap3$index_WAA_pointers = asap3$index_WAA_pointers[which_indices] - asap3$IAA_mats <- asap3$IAA_mats[which_indices] - asap3$use_survey <- asap3$use_index[which_indices] - data$n_indices <- asap3$n_indices + else { + asap3 = input$asap3 + which_indices <- which(asap3$use_index ==1) + asap3$n_indices = length(which_indices) + asap3$survey_index_units <- asap3$index_units[which_indices] + asap3$survey_acomp_units <- asap3$index_acomp_units[which_indices] + asap3$survey_WAA_pointers <- asap3$index_WAA_pointers[which_indices] + asap3$survey_month <- matrix(asap3$index_month[which_indices], asap3$n_years, asap3$n_indices, byrow = TRUE) + asap3$use_survey_acomp <- asap3$use_index_acomp[which_indices] + asap3$index_WAA_pointers = asap3$index_WAA_pointers[which_indices] + asap3$IAA_mats <- asap3$IAA_mats[which_indices] + asap3$use_survey <- asap3$use_index[which_indices] + data$n_indices <- asap3$n_indices } - data$agg_indices = matrix(NA, data$n_years_model, data$n_indices) - data$use_indices = matrix(1, data$n_years_model, data$n_indices) - data$agg_index_sigma = matrix(NA, data$n_years_model, data$n_indices) - data$index_paa = array(NA, dim = c(data$n_indices, data$n_years_model, data$n_ages)) - data$use_index_paa = matrix(1, data$n_years_model, data$n_indices) - data$index_Neff = matrix(NA, data$n_years_model, data$n_indices) + data$agg_indices = matrix(NA, data$n_years_model, data$n_indices) + data$use_indices = matrix(1, data$n_years_model, data$n_indices) + data$agg_index_sigma = matrix(NA, data$n_years_model, data$n_indices) + data$index_paa = array(NA, dim = c(data$n_indices, data$n_years_model, data$n_ages)) + data$use_index_paa = matrix(1, data$n_years_model, data$n_indices) + data$index_Neff = matrix(NA, data$n_years_model, data$n_indices) if(!is.null(asap3)) { data$units_indices <- asap3$survey_index_units data$fracyr_indices = (asap3$survey_month-1)/12 #make sure that this is right - for(i in 1:data$n_indices) data$agg_indices[,i] = asap3$IAA_mats[[i]][,2] for(i in 1:data$n_indices) { + data$agg_indices[,i] = asap3$IAA_mats[[i]][,2] for(y in 1:data$n_years_model) if(asap3$IAA_mats[[i]][y,2] < 1e-15) data$use_indices[y,i] = 0 } for(i in 1:data$n_indices) data$agg_index_sigma[,i] = asap3$IAA_mats[[i]][,3] @@ -43,15 +43,18 @@ set_indices = function(input, index_opts=NULL) temp[which(temp<0)] = 0 data$index_paa[i,,] = temp/apply(temp,1,sum) } + data$index_paa[is.na(data$index_paa)] = 0 for(i in 1:data$n_indices) { - if(asap3$use_survey_acomp[i] != 1) data$use_index_paa[,i] = 0 - else for(y in 1:data$n_years_model) if(asap3$IAA_mats[[i]][y,4 + data$n_ages] < 1e-15 | sum(data$index_paa[i,y,] > 1e-15)<2) data$use_index_paa[y,i] = 0 + if(asap3$use_survey_acomp[i] != 1){ + data$use_index_paa[,i] = 0 + } else { + for(y in 1:data$n_years_model) if(asap3$IAA_mats[[i]][y,4 + data$n_ages] < 1e-15 | sum(data$index_paa[i,y,] > 1e-15) < 2) data$use_index_paa[y,i] = 0 + } } data$units_index_paa <- asap3$survey_acomp_units for(i in 1:data$n_indices) data$index_Neff[,i] = asap3$IAA_mats[[i]][,4 + data$n_ages] - - data$selblock_pointer_indices = matrix(rep(asap3$n_fleet_sel_blocks + 1:data$n_indices, each = data$n_years_model), data$n_years_model, data$n_indices) + data$selblock_pointer_indices = matrix(rep(asap3$n_fleet_sel_blocks + 1:data$n_indices, each = data$n_years_model), data$n_years_model, data$n_indices) } else { diff --git a/R/set_selectivity.R b/R/set_selectivity.R index 47213cdb..8260acbd 100644 --- a/R/set_selectivity.R +++ b/R/set_selectivity.R @@ -32,8 +32,8 @@ set_selectivity = function(input, selectivity) } no_asap = is.null(asap3) selopts <- c("age-specific","logistic","double-logistic","decreasing-logistic") - if(!no_asap) data$n_selblocks <- asap3$n_fleet_sel_blocks + asap3$n_indices - else data$n_selblocks = data$n_fleets + data$n_indices + # if(!no_asap) data$n_selblocks <- asap3$n_fleet_sel_blocks + asap3$n_indices + # if(no_asap) data$n_selblocks = data$n_fleets + data$n_indices if(is.null(selectivity$model)) { #if(!no_asap) data$selblock_models <- c(asap3$sel_block_option, asap3$index_sel_option) @@ -145,9 +145,9 @@ set_selectivity = function(input, selectivity) # For age-specific selectivity blocks, check for ages with ~zero catch and fix these at 0 age_specific <- which(data$selblock_models==1) for(b in age_specific){ -# if(all(phase_selpars[b,] < 0)){ # if no selpars estimated, keep fixed at specified initial values -# phase_selpars[b,] = -1 -# } else { + if(all(phase_selpars[b,] < 0)){ # if no selpars estimated, keep fixed at specified initial values + phase_selpars[b,] = -1 + } else { ind = list(fleets = which(apply(data$selblock_pointer_fleets == b,2,sum) > 0)) ind$indices = which(apply(data$selblock_pointer_indices == b,2,sum) > 0) paa = matrix(nrow = 0, ncol = data$n_ages) @@ -164,7 +164,7 @@ set_selectivity = function(input, selectivity) y = apply(paa,2,sum) selpars_ini[b, which(y < 1e-5)] = 0 phase_selpars[b, which(y < 1e-5)] = -1 -# } + } } data$selpars_est <- phase_selpars data$selpars_est[data$selpars_est == -1] = 0