Skip to content

Commit

Permalink
fix errors in prepare_wham_input when missing indices or paa (butterf…
Browse files Browse the repository at this point in the history 
…ish, winter flounder)
  • Loading branch information
@brianstock-NOAA
brianstock-NOAA committed Aug 25, 2021
1 parent c782951 commit 12abdef
Show file tree
Hide file tree
Showing 3 changed files with 101 additions and 99 deletions.
133 changes: 66 additions & 67 deletions R/prepare_wham_input.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,115 +275,114 @@ 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
{
#if no asap3 is provided, make some default values to
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))
}
Expand Down
55 changes: 29 additions & 26 deletions R/set_indices.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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
{
Expand Down
12 changes: 6 additions & 6 deletions R/set_selectivity.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down

0 comments on commit 12abdef

Please sign in to comment.