updates

Author: Ovec8hkin

R/data_processing.R

@@ -401,8 +401,8 @@ get_reference_points <- function(model_runs, extra_columns, hcr_filter, om_filte
 get_b40_timeseries <- function(model_runs, extra_columns, hcr_filter, om_filter){
 
     get_rps <- function(om_name, hcr_name, recruitment, prop_fs){
-        om <- om_list[which(om_names == om_name)]
-        hcr <- hcr_list[which(hcr_names == hcr_name)]
+        om <- om_list[which(publication_oms == om_name)]
+        hcr <- hcr_list[which(publication_hcrs == hcr_name)]
 
         om <- om[[1]]
 
@@ -425,13 +425,13 @@ get_b40_timeseries <- function(model_runs, extra_columns, hcr_filter, om_filter)
         return(ref_pts$Bref)
     }
 
-    avg_recruitment <- get_recruits(model_runs, extra_columns2, hcr_filter, om_filter) %>% 
+    avg_recruitment <- get_recruits(model_runs, extra_columns, hcr_filter, om_filter) %>% 
         filter(L1 == "naa") %>% 
         group_by(sim, om, hcr) %>%
         mutate(avg_rec = unlist(lapply(slide(rec, ~.x, .before=Inf), \(x) mean(x)))) %>%
         arrange(hcr, om, sim)
 
-    prop_fs_df <- get_fishing_mortalities(model_runs, extra_columns2, hcr_filter, om_filter) %>%
+    prop_fs_df <- get_fishing_mortalities(model_runs, extra_columns, hcr_filter, om_filter) %>%
         filter(L1 != "faa_est") %>%
         group_by(time, sim, om, hcr, fleet) %>%
         mutate(

R/data_utils.R

@@ -72,13 +72,16 @@ relativize_performance <- function(data, rel_column, value_column, rel_value, gr
     if(is.null(rel_value)){
         return(data)
     }
+
+    total_cols <- ncol(data)-1
 
     return(
         data %>%
             group_by(across(all_of(grouping))) %>%
             pivot_wider(names_from=rel_column, values_from = value_column) %>%
             mutate(across(everything(), ~ . / eval(rlang::parse_expr(rel_value)))) %>%
-            pivot_longer((length(grouping)+1):(ncol(.)), names_to=rel_column, values_to=value_column)
+            # pivot_longer((length(grouping)+1):(ncol(.)), names_to=rel_column, values_to=value_column)
+            pivot_longer(total_cols:(ncol(.)), names_to=rel_column, values_to=value_column)
     )
 }
 
@@ -175,7 +178,7 @@ scale_and_rank <- function(data, col_name){
                 rank = ifelse(
                             name %in% c(
                                 "Catch AAV", 
-                                "Proportion of Years with Low SSB", 
+                                "Proportion of Years SSB < B35", 
                                 "Recovery Time"
                             ), 
                             factor(desc(row_number())), 

R/performance_metrics.R

@@ -321,11 +321,12 @@ prop_low_biomass <- function(
             total_low_bio = sum(low_bio),
             prop_years = total_low_bio/(time_horizon[2]-time_horizon[1])
         ) %>%
+        select(-c("total_low_bio")) %>%
         relativize_performance(
             rel_column = "hcr",
             value_column = "prop_years",
             rel_value = relative,
-            grouping = group_columns
+            grouping = c("sim", "om")
         )
 
     if(!is.null(extra_filter)){
@@ -915,7 +916,8 @@ average_annual_value <- function(
     time_horizon = c(65, NA), 
     extra_filter=NULL, 
     relative=NULL, 
-    summarise_by=c("om", "hcr")
+    summarise_by=c("om", "hcr"),
+    summary_out=TRUE
 ){
 
     group_columns <- c("sim", summarise_by)
@@ -953,11 +955,14 @@ average_annual_value <- function(
         avg_rel_value <- avg_rel_value %>% filter(eval(extra_filter))
     }
 
-    return(
-        avg_rel_value %>%
+    out <- avg_rel_value
+    if(summary_out){
+        out <- avg_rel_value %>%
             group_by(across(all_of(summarise_by))) %>%
             median_qi(annual_value, .width=interval_widths, .simple_names=FALSE)
-    )
+    }
+
+    return(out)
 }
 
 #' Compute Average Annual Dynamic Value of Catch across projection period
@@ -992,7 +997,8 @@ average_annual_dynamic_value <- function(
     time_horizon = c(65, NA), 
     extra_filter=NULL, 
     relative=NULL, 
-    summarise_by=c("om", "hcr")
+    summarise_by=c("om", "hcr"),
+    summary_out=TRUE
 ){
 
     compute_dynamic_value <- function(landings, min_price_age, max_price_age, breakpoints=c(15, 30)){
@@ -1057,11 +1063,14 @@ average_annual_dynamic_value <- function(
         dyn_value <- dyn_value %>% filter(eval(extra_filter))
     }
 
-    return(
-        dyn_value %>%
+    out <- dyn_value
+    if(summary_out){
+        out <- dyn_value %>%
             group_by(across(all_of(summarise_by))) %>%
             median_qi(dyn_annual_value, .width=interval_widths, .simple_names=FALSE)
-    )
+    }
+
+    return(out)
 
 }
 

R/plotting.R

@@ -319,17 +319,17 @@ plot_ssb_catch <- function(ssb_data, catch_data, v1="hcr", v2=NA, v3=NA, common_
 
     plot <- ggplot(d) + 
         geom_line(data = base_hcr, aes(x=time, y=median, group=.data[[v1]], color=.data[[v1]]), size=0.85)+
-        geom_line(aes(x=time, y=median, group=.data[[v1]], color=color_group, size=color_group))+
-        geom_line(
-            data = d %>% filter(color_group != "Other", time > common_trajectory-1), 
-            aes(x=time, y=median, ymin=lower, ymax=upper, group=.data[[v1]], color=color_group, size=color_group)
-        )+
+        geom_line(aes(x=time, y=median, group=.data[[v1]], color=hcr), size=0.85)+
+        # geom_line(
+        #     data = d %>% filter(color_group != "Other", time > common_trajectory-1), 
+        #     aes(x=time, y=median, ymin=lower, ymax=upper, group=.data[[v1]], color=color_group, size=color_group)
+        # )+
         geom_line(data = common, aes(x=time, y=median), size=0.85)+
         geom_vline(data=common, aes(xintercept=common), linetype="dashed")+
         # geom_hline(yintercept=121.4611, linetype="dashed")+
         scale_fill_brewer(palette="Blues")+
         scale_color_manual(values=colors)+
-        scale_size_manual(values=sizes)+
+        # scale_size_manual(values=sizes)+
         # scale_y_continuous(limits=c(0, 320))+
         labs(x="Year", y="1000s Metric Tons")+
         coord_cartesian(expand=0)+
@@ -666,7 +666,23 @@ plot_performance_metric_summary <- function(perf_data, v1="hcr", v2="om", is_rel
     metric_minmax = perf_data %>% group_by(name) %>% summarise(min=min(lower), max=max(upper))
     axis_scalar <- c(0.9, 1.1)
 
+    # hist_abcs <- c(44200, 37100, 33400, 28800, 25200, 25000, 28800, 25300, 19600, 17200, 16800, 15900, 17200, 16900, 17300, 20900, 23000, 21000, 21000, 20100, 18000, 16100, 15200, 16000, 17200, 16200, 13700, 13700, 11800, 13100, 15000, 15100, 22000, 29600, 34500, 40500)
+
     summary <- perf_data %>% filter(!is.infinite(median), hcr != "No Fishing") %>% summarise(median=mean(median))
+    # summary$median <- rep(
+    #     c(
+    #         18538/1000,#median(assessment$t.series[,"Catch_HAL"]+assessment$t.series[,"Catch_TWL"]),
+    #         perf_data %>% ungroup() %>% filter(!is.infinite(median), hcr != "No Fishing", name=="Catch AAV") %>% summarise(median=mean(median)) %>% as.numeric,
+    #         # median(assessment$t.series[,"spbiom"]),
+    #         105.935,
+    #         median(apply(assessment$natage.female, 1, \(x) compute_average_age(x, 2:31))),
+    #         sum(assessment$t.series[,"spbiom"] < 105)/length(assessment$t.series[,"spbiom"])
+    #     ),
+    #     length(perf_data$om %>% unique)*2
+    # )
+    if(is_relative){
+        summary$median <- rep(1, nrow(summary))
+    }
 
     perf_data <- perf_data %>% mutate(color_group = as.character(hcr))
     if(!is.null(highlight)){
@@ -691,12 +707,12 @@ plot_performance_metric_summary <- function(perf_data, v1="hcr", v2="om", is_rel
     color_var <- ifelse(is.null(highlight), "rank", "color_group")
 
     plot <- ggplot(perf_data)+
-                geom_vline(data=summary, aes(xintercept = median), color="black")+
+                geom_vline(data=summary, aes(xintercept = median), color="grey50", linetype="dashed")+
                 scale_shape_discrete()+
                 scale_color_manual(values=colors)+
                 # scale_color_manual(values=hcr_colors)+
                 # facet_wrap(vars(name), scales="free_x")+
-                labs(y="", x="", shape="OM", color="Performance Order")+
+                labs(y="", x="", shape="OM", color="Relative MS Order")+
                 coord_cartesian(expand=0)+
                 guides(shape="none", color=guide_legend(nrow=1))+
                 theme(
@@ -721,7 +737,7 @@ plot_performance_metric_summary <- function(perf_data, v1="hcr", v2="om", is_rel
                 ggh4x::facetted_pos_scales(
                     x = list(
                         scale_x_continuous(limits=c(0, 55)),
-                        scale_x_continuous(limits=c(0, 0.06), breaks=c(0, 0.02, 0.04, 0.06)),
+                        scale_x_continuous(limits=c(0, 0.08), breaks=c(0, 0.025, 0.05, 0.075)),
                         # scale_x_continuous(limits=c(0, 1), breaks=c(0, 0.50, 1.0)),
                         scale_x_continuous(limits=c(0, 550), breaks=c(0, 150, 300, 450)),
                         scale_x_continuous(limits=c(0, 15)),

R/run_mse.R

@@ -141,7 +141,7 @@ run_mse <- function(om, mp, mse_options, nyears_input=NA, seed=1120, file_suffix
         }
 
         prev_naa <- naa[y,,,, drop = FALSE]
-        out_vars <- project(
+        out_vars <- project_single(
             removals = removals_input,
             dem_params=dp_y,
             prev_naa=prev_naa,

R/run_mse_multiple.R

@@ -50,7 +50,7 @@ run_mse_multiple <- function(om_list, hcr_list, seed_list, mse_options_list, nye
             mse_run <- run_mse_parallel(nsims2, seeds, om, hcr, mse_options=opt, nyears=nyears, diagnostics=diagnostics)
 
             counter <- counter+1
-            mse_objects[[j]] <- mse_run
+            mse_objects[[ifelse(nsim_iters > 1, j, i)]] <- mse_run
 
             print(save)
             if(save || nsim_iters>1){
@@ -72,4 +72,4 @@ run_mse_multiple <- function(om_list, hcr_list, seed_list, mse_options_list, nye
 
     return(mse_objects)
 
-}
+}