diff --git a/R/RatioBiasEM.R b/R/RatioBiasEM.R index 689ebb1..2fb179d 100644 --- a/R/RatioBiasEM.R +++ b/R/RatioBiasEM.R @@ -15,7 +15,7 @@ get_RatioEM_catch_df<-function(EM_dir, dat, dat_yrs, #EM_dir = EM_out_dir; dat = new_EM_dat changeup=10, changedown=10,# add inputs, min & max allowable annual percent change in catch - c1=0.75,# add inputs, multiplicative constant + c1=1,# add inputs, multiplicative constant option="default", ...) { # B_ratio_denominator @@ -308,8 +308,8 @@ get_RatioEM_catch_df<-function(EM_dir, dat, dat_yrs, RatioBiasEM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = FALSE, - nyrs_assess, dat_yrs, sample_struct = NULL, sample_struct_hist = NULL, - seed = NULL, OM_out_dir, ...) { + nyrs_assess, dat_yrs, sample_struct = NULL, sample_struct_hist = NULL, + seed = NULL, OM_out_dir, ...) { SSMSE:::check_dir(EM_out_dir) # TODO: change this name to make it less ambiguous new_datfile_name <- "init_dat.ss" @@ -489,7 +489,7 @@ RatioBiasEM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = F new_OM_catch_list$discards$Discard <- new_OM_catch_list$discards$Discard * tmp_discards$bias } if(2 %in% OM_dat$fleetinfo$type){ # if bycatch fleet -- remove from catch list and keep only bycatch fleets in catch_F list - + byc_f <- which(OM_dat$fleetinfo$type==2)#as.numeric(row.names(OM_dat$fleetinfo[which(OM_dat$fleetinfo$type==2),])) new_OM_catch_list$catch<- new_OM_catch_list$catch[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),] new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),] @@ -497,7 +497,7 @@ RatioBiasEM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = F # new_OM_catch_list$catch<- new_OM_catch_list$catch[new_OM_catch_list$catch$fleet!=byc_f,] # new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[new_OM_catch_list$catch$fleet!=byc_f,] # new_OM_catch_list$catch_F<- new_OM_catch_list$catch_F[new_OM_catch_list$catch_F$fleet==byc_f,] - + } else{ new_OM_catch_list$catch_F <- NULL } @@ -508,7 +508,408 @@ RatioBiasEM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = F } - +RatioBiasEM2 <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = FALSE, + nyrs_assess, dat_yrs, sample_struct = NULL, sample_struct_hist = NULL, + seed = NULL, OM_out_dir, ...) { + SSMSE:::check_dir(EM_out_dir) + # TODO: change this name to make it less ambiguous + new_datfile_name <- "init_dat.ss" + # change the name of data file. + start <- SS_readstarter(file.path(EM_out_dir, "starter.ss"), + verbose = FALSE + ) + + if (init_loop) { + + # copy over raw data file from the OM to EM folder + SS_writedat(OM_dat, + file.path(EM_out_dir, new_datfile_name), + overwrite = TRUE, + verbose = FALSE + ) + orig_datfile_name <- start[["datfile"]] # save the original data file name. + start[["datfile"]] <- new_datfile_name + start[["seed"]] <- seed + SS_writestarter(start, file.path(EM_out_dir), + verbose = FALSE, + overwrite = TRUE + ) + ### NOTE:: BY THIS POINT, THE EM DATAFILE HAS BEEN UPDATED TO INCLUDE FORECAST YRS. + + # make sure the data file has the correct formatting (use existing data + # file in the EM directory to make sure)?? + # TODO: is this necessary, given we have sample structures? + new_EM_dat <- biasEM_change_dat( # edited func to account for historical bias + OM_datfile = new_datfile_name, + EM_datfile = orig_datfile_name, + EM_dir = EM_out_dir, + do_checks = TRUE, + verbose = verbose, + sample_struct_hist = sample_struct_hist + ) + ctl <- SS_readctl(file.path(EM_out_dir, start[["ctlfile"]]), + datlist = new_EM_dat + ) + if (ctl[["EmpiricalWAA"]] == 1) { + message( + "EM uses weight at age, so copying over wtatage file from OM.", + "\nNote wtatage data is not sampled." + ) + file.copy( + from = file.path(OM_out_dir, "wtatage.ss"), + to = file.path(EM_out_dir, "wtatage.ss"), + overwrite = TRUE + ) + } + if (!all(ctl[["time_vary_auto_generation"]] == 1)) { + warning("Turning off autogeneration of time varying lines in the control file of the EM") + ctl[["time_vary_auto_generation"]] <- rep(1, times = 5) + r4ss::SS_writectl(ctl, file.path(EM_out_dir, start[["ctlfile"]]), + overwrite = TRUE + ) + } + + } else { # end if init_loop==T + if (!is.null(sample_struct)) { + sample_struct_sub <- lapply(sample_struct, + function(df, y) df[df[, 1] %in% y, ], + y = dat_yrs - nyrs_assess + ) + } else { + sample_struct_sub <- NULL + } + + new_EM_dat <- add_new_dat_BIAS( ######## NEW FUNCTION TO BUILD IN CONVERSION FACTOR + OM_dat = OM_dat, + EM_datfile = new_datfile_name, + sample_struct = sample_struct_sub, + EM_dir = EM_out_dir, + nyrs_assess = nyrs_assess, + do_checks = TRUE, + new_datfile_name = new_datfile_name, + verbose = verbose + ) + + } # end else not first iteration + + # Update SS random seed + start <- SS_readstarter(file.path(EM_out_dir, "starter.ss"), + verbose = FALSE + ) + start[["seed"]] <- seed + SS_writestarter(start, file.path(EM_out_dir), + verbose = FALSE, + overwrite = TRUE + ) + # manipulate the forecasting file. + # make sure enough yrs can be forecasted. + + fcast <- SS_readforecast(file.path(EM_out_dir, "forecast.ss"), + readAll = TRUE, + verbose = FALSE + ) + # check that it can be used in the EM. fleets shoul + SSMSE:::check_EM_forecast(fcast, + n_flts_catch = length(which(new_EM_dat[["fleetinfo"]][, "type"] %in% + c(1, 2))) + ) + fcast <- SSMSE:::change_yrs_fcast(fcast, + make_yrs_rel = (init_loop == TRUE), + nyrs_fore = nyrs_assess, + mod_styr = new_EM_dat[["styr"]], + mod_endyr = new_EM_dat[["endyr"]] + ) + SS_writeforecast(fcast, + dir = EM_out_dir, writeAll = TRUE, overwrite = TRUE, + verbose = FALSE + ) + # given all checks are good, run the EM + # check convergence (figure out way to error if need convergence) + # get the future catch using the management strategy used in the SS model. + run_EM(EM_dir = EM_out_dir, verbose = verbose, check_converged = TRUE) + # get the forecasted catch. + new_EM_catch_list <- get_RatioEM_catch_df(EM_dir = EM_out_dir, dat = new_EM_dat, dat_yrs=dat_yrs, + changeup=10, changedown=10,c1=0.75) + ## COnSIDER MAKING A get_OM_catch_df if structure of OM =/= EM. + # For the simple approach, we can just apply a series of scalars from the EM catch list to create an OM catch list + new_OM_catch_list = new_EM_catch_list + + ## IF FixedCatches==TRUE + # Will need to be mindful about units -- so far, assuming is presented in same units as historical OM + # come back to this section if want to use different units + if(!is.null(sample_struct$FixedCatch)){ + + # tmp_ss<- sample_struct$FixedCatch[sample_struct$FixedCatch$year==dat_yrs,] + tmp_ss<- sample_struct$FixedCatch[sample_struct$FixedCatch$year %in% dat_yrs,] # create obj of fixed catches + colnames(tmp_ss)[4]<- "Fcatch" # rename fixed catches to allow for merge + + if(nrow(tmp_ss)>0){ + if(!is.null(new_OM_catch_list$catch)){ + tmp_ss_catch <- tmp_ss[tmp_ss$units!=99,] + if(nrow(tmp_ss_catch)>0){ + tmp_merge <- base::merge(base::abs(new_OM_catch_list$catch), base::abs(tmp_ss_catch), all.x=TRUE, all.y=FALSE) # merge + tmp_merge$catch[which(!is.na(tmp_merge$Fcatch))] <- tmp_merge$Fcatch[which(!is.na(tmp_merge$Fcatch))] # replace fixed catches with + new_OM_catch_list$catch<-tmp_merge[,c(1:5)] #reorder columns of merged + } + }#end if catch exists + + if(!is.null(new_OM_catch_list$catch_bio)){ + tmp_ss_catch <- tmp_ss[tmp_ss$units==1,] + if(nrow(tmp_ss_catch)>0){ + tmp_merge <- base::merge(base::abs(new_OM_catch_list$catch), base::abs(tmp_ss_catch), all.x=TRUE, all.y=FALSE) # merge + tmp_merge$catch[which(!is.na(tmp_merge$Fcatch))] <- tmp_merge$Fcatch[which(!is.na(tmp_merge$Fcatch))] # replace fixed catches with + new_OM_catch_list$catch_bio<-tmp_merge[,c(1:5)] #reorder columns of merged + } + }else{ + new_OM_catch_list$catch_bio <- NULL }#end if catch_bio exists + + if(!is.null(new_OM_catch_list$catch_F)){ + tmp_ss_catch <- tmp_ss[tmp_ss$units==99,] + if(nrow(tmp_ss_catch)>0){ + tmp_merge <- base::merge(base::abs(new_OM_catch_list$catch_F), base::abs(tmp_ss_catch), all.x=TRUE, all.y=FALSE) # merge + tmp_merge$catch[which(!is.na(tmp_merge$Fcatch))] <- tmp_merge$Fcatch[which(!is.na(tmp_merge$Fcatch))] # replace fixed catches with + new_OM_catch_list$catch_F<-tmp_merge[,c(1:5)] #reorder columns of merged + } + }#end if catch exists + }# end if fixed catches in this mgmt cycle. + + }# end fixed catches + + ## EM2OM bias correction + if(!is.null(new_OM_catch_list$catch)){ + tmp_catch <- base::merge(new_OM_catch_list$catch, sample_struct$EM2OMcatch_bias, all.x=TRUE, all.y=FALSE) + new_OM_catch_list$catch$catch <- new_OM_catch_list$catch$catch * tmp_catch$bias + } + if(!is.null(new_OM_catch_list$catch_bio)){ + tmp_catch_bio <- base::merge(new_OM_catch_list$catch_bio, sample_struct$EM2OMcatch_bias, all.x=TRUE) + new_OM_catch_list$catch_bio$catch <- new_OM_catch_list$catch_bio$catch * tmp_catch_bio$bias + } + # new_OM_catch_list$catch_bio <- NULL + if(!is.null(new_OM_catch_list$discards)){ + tmp_discards <- base::merge(base::abs(new_OM_catch_list$discards), base::abs(sample_struct$EM2OMdiscard_bias), all.x=TRUE) + tmp_discards <- tmp_discards[base::order(base::abs(tmp_discards$Flt),base::abs(tmp_discards$Yr),base::abs(tmp_discards$Seas)),] + new_OM_catch_list$discards$Discard <- new_OM_catch_list$discards$Discard * tmp_discards$bias + } + if(2 %in% OM_dat$fleetinfo$type){ # if bycatch fleet -- remove from catch list and keep only bycatch fleets in catch_F list + + byc_f <- which(OM_dat$fleetinfo$type==2)#as.numeric(row.names(OM_dat$fleetinfo[which(OM_dat$fleetinfo$type==2),])) + new_OM_catch_list$catch<- new_OM_catch_list$catch[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),] + new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),] + new_OM_catch_list$catch_F<- new_OM_catch_list$catch_F[which(is.element(new_OM_catch_list$catch_F$fleet,byc_f)),] + # new_OM_catch_list$catch<- new_OM_catch_list$catch[new_OM_catch_list$catch$fleet!=byc_f,] + # new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[new_OM_catch_list$catch$fleet!=byc_f,] + # new_OM_catch_list$catch_F<- new_OM_catch_list$catch_F[new_OM_catch_list$catch_F$fleet==byc_f,] + + } else{ + new_OM_catch_list$catch_F <- NULL + } + + + new_catch_list<-new_OM_catch_list + new_catch_list + +} +RatioBiasEM3 <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = FALSE, + nyrs_assess, dat_yrs, sample_struct = NULL, sample_struct_hist = NULL, + seed = NULL, OM_out_dir, ...) { + SSMSE:::check_dir(EM_out_dir) + # TODO: change this name to make it less ambiguous + new_datfile_name <- "init_dat.ss" + # change the name of data file. + start <- SS_readstarter(file.path(EM_out_dir, "starter.ss"), + verbose = FALSE + ) + + if (init_loop) { + + # copy over raw data file from the OM to EM folder + SS_writedat(OM_dat, + file.path(EM_out_dir, new_datfile_name), + overwrite = TRUE, + verbose = FALSE + ) + orig_datfile_name <- start[["datfile"]] # save the original data file name. + start[["datfile"]] <- new_datfile_name + start[["seed"]] <- seed + SS_writestarter(start, file.path(EM_out_dir), + verbose = FALSE, + overwrite = TRUE + ) + ### NOTE:: BY THIS POINT, THE EM DATAFILE HAS BEEN UPDATED TO INCLUDE FORECAST YRS. + + # make sure the data file has the correct formatting (use existing data + # file in the EM directory to make sure)?? + # TODO: is this necessary, given we have sample structures? + new_EM_dat <- biasEM_change_dat( # edited func to account for historical bias + OM_datfile = new_datfile_name, + EM_datfile = orig_datfile_name, + EM_dir = EM_out_dir, + do_checks = TRUE, + verbose = verbose, + sample_struct_hist = sample_struct_hist + ) + ctl <- SS_readctl(file.path(EM_out_dir, start[["ctlfile"]]), + datlist = new_EM_dat + ) + if (ctl[["EmpiricalWAA"]] == 1) { + message( + "EM uses weight at age, so copying over wtatage file from OM.", + "\nNote wtatage data is not sampled." + ) + file.copy( + from = file.path(OM_out_dir, "wtatage.ss"), + to = file.path(EM_out_dir, "wtatage.ss"), + overwrite = TRUE + ) + } + if (!all(ctl[["time_vary_auto_generation"]] == 1)) { + warning("Turning off autogeneration of time varying lines in the control file of the EM") + ctl[["time_vary_auto_generation"]] <- rep(1, times = 5) + r4ss::SS_writectl(ctl, file.path(EM_out_dir, start[["ctlfile"]]), + overwrite = TRUE + ) + } + + } else { # end if init_loop==T + if (!is.null(sample_struct)) { + sample_struct_sub <- lapply(sample_struct, + function(df, y) df[df[, 1] %in% y, ], + y = dat_yrs - nyrs_assess + ) + } else { + sample_struct_sub <- NULL + } + + new_EM_dat <- add_new_dat_BIAS( ######## NEW FUNCTION TO BUILD IN CONVERSION FACTOR + OM_dat = OM_dat, + EM_datfile = new_datfile_name, + sample_struct = sample_struct_sub, + EM_dir = EM_out_dir, + nyrs_assess = nyrs_assess, + do_checks = TRUE, + new_datfile_name = new_datfile_name, + verbose = verbose + ) + + } # end else not first iteration + + # Update SS random seed + start <- SS_readstarter(file.path(EM_out_dir, "starter.ss"), + verbose = FALSE + ) + start[["seed"]] <- seed + SS_writestarter(start, file.path(EM_out_dir), + verbose = FALSE, + overwrite = TRUE + ) + # manipulate the forecasting file. + # make sure enough yrs can be forecasted. + + fcast <- SS_readforecast(file.path(EM_out_dir, "forecast.ss"), + readAll = TRUE, + verbose = FALSE + ) + # check that it can be used in the EM. fleets shoul + SSMSE:::check_EM_forecast(fcast, + n_flts_catch = length(which(new_EM_dat[["fleetinfo"]][, "type"] %in% + c(1, 2))) + ) + fcast <- SSMSE:::change_yrs_fcast(fcast, + make_yrs_rel = (init_loop == TRUE), + nyrs_fore = nyrs_assess, + mod_styr = new_EM_dat[["styr"]], + mod_endyr = new_EM_dat[["endyr"]] + ) + SS_writeforecast(fcast, + dir = EM_out_dir, writeAll = TRUE, overwrite = TRUE, + verbose = FALSE + ) + # given all checks are good, run the EM + # check convergence (figure out way to error if need convergence) + # get the future catch using the management strategy used in the SS model. + run_EM(EM_dir = EM_out_dir, verbose = verbose, check_converged = TRUE) + # get the forecasted catch. + new_EM_catch_list <- get_RatioEM_catch_df(EM_dir = EM_out_dir, dat = new_EM_dat, dat_yrs=dat_yrs, + changeup=0.2, changedown=0.2,c1=1) + ## COnSIDER MAKING A get_OM_catch_df if structure of OM =/= EM. + # For the simple approach, we can just apply a series of scalars from the EM catch list to create an OM catch list + new_OM_catch_list = new_EM_catch_list + + ## IF FixedCatches==TRUE + # Will need to be mindful about units -- so far, assuming is presented in same units as historical OM + # come back to this section if want to use different units + if(!is.null(sample_struct$FixedCatch)){ + + # tmp_ss<- sample_struct$FixedCatch[sample_struct$FixedCatch$year==dat_yrs,] + tmp_ss<- sample_struct$FixedCatch[sample_struct$FixedCatch$year %in% dat_yrs,] # create obj of fixed catches + colnames(tmp_ss)[4]<- "Fcatch" # rename fixed catches to allow for merge + + if(nrow(tmp_ss)>0){ + if(!is.null(new_OM_catch_list$catch)){ + tmp_ss_catch <- tmp_ss[tmp_ss$units!=99,] + if(nrow(tmp_ss_catch)>0){ + tmp_merge <- base::merge(base::abs(new_OM_catch_list$catch), base::abs(tmp_ss_catch), all.x=TRUE, all.y=FALSE) # merge + tmp_merge$catch[which(!is.na(tmp_merge$Fcatch))] <- tmp_merge$Fcatch[which(!is.na(tmp_merge$Fcatch))] # replace fixed catches with + new_OM_catch_list$catch<-tmp_merge[,c(1:5)] #reorder columns of merged + } + }#end if catch exists + + if(!is.null(new_OM_catch_list$catch_bio)){ + tmp_ss_catch <- tmp_ss[tmp_ss$units==1,] + if(nrow(tmp_ss_catch)>0){ + tmp_merge <- base::merge(base::abs(new_OM_catch_list$catch), base::abs(tmp_ss_catch), all.x=TRUE, all.y=FALSE) # merge + tmp_merge$catch[which(!is.na(tmp_merge$Fcatch))] <- tmp_merge$Fcatch[which(!is.na(tmp_merge$Fcatch))] # replace fixed catches with + new_OM_catch_list$catch_bio<-tmp_merge[,c(1:5)] #reorder columns of merged + } + }else{ + new_OM_catch_list$catch_bio <- NULL }#end if catch_bio exists + + if(!is.null(new_OM_catch_list$catch_F)){ + tmp_ss_catch <- tmp_ss[tmp_ss$units==99,] + if(nrow(tmp_ss_catch)>0){ + tmp_merge <- base::merge(base::abs(new_OM_catch_list$catch_F), base::abs(tmp_ss_catch), all.x=TRUE, all.y=FALSE) # merge + tmp_merge$catch[which(!is.na(tmp_merge$Fcatch))] <- tmp_merge$Fcatch[which(!is.na(tmp_merge$Fcatch))] # replace fixed catches with + new_OM_catch_list$catch_F<-tmp_merge[,c(1:5)] #reorder columns of merged + } + }#end if catch exists + }# end if fixed catches in this mgmt cycle. + + }# end fixed catches + + ## EM2OM bias correction + if(!is.null(new_OM_catch_list$catch)){ + tmp_catch <- base::merge(new_OM_catch_list$catch, sample_struct$EM2OMcatch_bias, all.x=TRUE, all.y=FALSE) + new_OM_catch_list$catch$catch <- new_OM_catch_list$catch$catch * tmp_catch$bias + } + if(!is.null(new_OM_catch_list$catch_bio)){ + tmp_catch_bio <- base::merge(new_OM_catch_list$catch_bio, sample_struct$EM2OMcatch_bias, all.x=TRUE) + new_OM_catch_list$catch_bio$catch <- new_OM_catch_list$catch_bio$catch * tmp_catch_bio$bias + } + # new_OM_catch_list$catch_bio <- NULL + if(!is.null(new_OM_catch_list$discards)){ + tmp_discards <- base::merge(base::abs(new_OM_catch_list$discards), base::abs(sample_struct$EM2OMdiscard_bias), all.x=TRUE) + tmp_discards <- tmp_discards[base::order(base::abs(tmp_discards$Flt),base::abs(tmp_discards$Yr),base::abs(tmp_discards$Seas)),] + new_OM_catch_list$discards$Discard <- new_OM_catch_list$discards$Discard * tmp_discards$bias + } + if(2 %in% OM_dat$fleetinfo$type){ # if bycatch fleet -- remove from catch list and keep only bycatch fleets in catch_F list + + byc_f <- which(OM_dat$fleetinfo$type==2)#as.numeric(row.names(OM_dat$fleetinfo[which(OM_dat$fleetinfo$type==2),])) + new_OM_catch_list$catch<- new_OM_catch_list$catch[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),] + new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),] + new_OM_catch_list$catch_F<- new_OM_catch_list$catch_F[which(is.element(new_OM_catch_list$catch_F$fleet,byc_f)),] + # new_OM_catch_list$catch<- new_OM_catch_list$catch[new_OM_catch_list$catch$fleet!=byc_f,] + # new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[new_OM_catch_list$catch$fleet!=byc_f,] + # new_OM_catch_list$catch_F<- new_OM_catch_list$catch_F[new_OM_catch_list$catch_F$fleet==byc_f,] + + } else{ + new_OM_catch_list$catch_F <- NULL + } + + + new_catch_list<-new_OM_catch_list + new_catch_list + +} #### add_new_dat Function Line-By-Line #### add_new_dat_BIAS<- function (OM_dat, EM_datfile, sample_struct, EM_dir, nyrs_assess,