Skip to content

Commit

Permalink
docs: update readme to use develop_OMs
Browse files Browse the repository at this point in the history 
For some reason, runs look really different than in the past, so will need to figure out why.
  • Loading branch information
@k-doering-NOAA
k-doering-NOAA committed Sep 14, 2020
1 parent 309260e commit ee336f5
Show file tree
Hide file tree
Showing 4 changed files with 29 additions and 126 deletions.
61 changes: 12 additions & 49 deletions README.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,47 +101,17 @@ dir.create(run_SSMSE_dir)
```
### Create the operating models (OMs)

The cod model with H = 0.65 (as in scenario 1) is included as external package data in SSMSE. However, we will need to modify it to use as an operating model with H = 1 (as in scenario 2).
The cod model with H = 0.65 (as in scenario 1) is included as external package data in SSMSE. However, we will need to modify it to use as an operating model with H = 1 (as in scenario 2). Note in this case that refit_OM is false, so the model is not being refit, just run through without fitting. To condition the new model on the same data as the input model, refit_OM should be TRUE.

```{r}
cod_mod_path <- system.file("extdata", "models", "cod", package = "SSMSE")
# copy to a new location:
file.copy(from = cod_mod_path, to = run_SSMSE_dir, recursive = TRUE)
file.rename(from = file.path(run_SSMSE_dir, "cod"), to = file.path(run_SSMSE_dir, "cod-1"))
cod_1_path <- file.path(run_SSMSE_dir, "cod-1")
# make model read initial values from control file and not ss.par
start <- r4ss::SS_readstarter(file = file.path(cod_1_path, "starter.ss"), verbose = FALSE)
start$init_values_src # verify reading from the control file
# change the steepness from 0.2 to 0.1 in the control files
r4ss::SS_changepars(dir = cod_1_path, ctlfile = "control.ss_new",
newctlfile = "control_modified.ss", strings = "SR_BH_steep", newvals = 1)
# read in the par file to save the historic recruitment deviations:
parfile <- r4ss::SS_readpar_3.30(parfile = file.path(cod_1_path, "ss.par"),
datsource = file.path(cod_1_path, "ss3.dat"),
ctlsource = file.path(cod_1_path, "control.ss_new"),
verbose = FALSE)
# remove files with old steepness values
file.remove(file.path(cod_1_path, "control.ss_new"))
file.remove(file.path(cod_1_path, "control.ss"))
file.remove(file.path(cod_1_path, "ss.par"))
file.rename(from = file.path(cod_1_path, "control_modified.ss"),
to = file.path(cod_1_path, "control.ss"))
```
Rerun this model with no estimation to get valid ss.par and control.ss_new files,
then add in the historical recruitment deviations:
```{r, results = "hide"}
SSMSE:::run_ss_model(dir = cod_1_path,
admb_options = "-maxfn 0 -phase 50 -nohess",
verbose = FALSE)
# add back original recdevs into the model (b/c not specified through the ctl file)
new_parfile <- r4ss::SS_readpar_3.30(parfile = file.path(cod_1_path, "ss.par"),
datsource = file.path(cod_1_path, "ss3.dat"),
ctlsource = file.path(cod_1_path, "control.ss"), verbose = FALSE)
# add in the recdevs to new the parfile
new_parfile$recdev1[, "recdev"] <- parfile$recdev1[, "recdev"]
r4ss::SS_writepar_3.30(new_parfile, outfile = file.path(cod_1_path, "ss.par"),
verbose = FALSE)
develop_OMs(OM_name = "cod", out_dir = run_SSMSE_dir, par_name = "SR_BH_steep",
par_vals = 1, refit_OMs = FALSE, hess = FALSE)
# OM model for scenario 2
cod_1_path <- file.path(run_SSMSE_dir, "cod_SR_BH_steep_1")
```

### Examine the management procedure used

We will use the same management procedure for both scenarios:
Expand Down Expand Up @@ -229,7 +199,7 @@ run_SSMSE(scen_name_vec = c("H-ctl", "H-1"),# name of the scenario
scope = "2", # to use the same recruitment devs across scenarios.
impl_error_pattern = "none", # Don't use implementation error
run_EM_last_yr = FALSE, # Run the EM in 106
run_parallel = TRUE,
run_parallel = FALSE,
sample_struct_list = sample_struct_list, # How to sample data for running the EM.
seed = 12345) #Set a fixed integer seed that allows replication
```
Expand All @@ -251,26 +221,19 @@ library(tidyr) # use install.packages("tidyr") to install package if needed
library(dplyr)
```
```{r plot_SSB}
summary$ts <- tidyr::separate(summary$ts,
col = model_run,
into = c(NA, "model_type"),
remove = FALSE,
sep = "_",
extra = "drop")
# check values for cod_OM
summary$scalar %>%
dplyr::filter(iteration == 1) %>%
dplyr::filter(scenario == "H-1") %>%
dplyr::select(iteration, scenario, SR_BH_steep, model_run)
# plot SSB by year and model run - need to correct using code from the
# think tank
ggplot2::ggplot(data = subset(summary$ts, model_run %in% c("cod_OM", "cod-1_OM", "cod_EM_106")),
# plot SSB by year and model run
ggplot2::ggplot(data = subset(summary$ts, model_run %in% c("cod_OM", "cod_SR_BH_steep_1_OM", "cod_EM_103")),
ggplot2::aes(x = year, y = SpawnBio)) +
ggplot2::geom_vline(xintercept = 100, color = "gray") +
ggplot2::geom_line(ggplot2::aes(linetype = as.character(iteration), color = model_type))+
ggplot2::scale_color_manual(values = c("#D65F00", "black")) +
ggplot2::geom_line(ggplot2::aes(linetype = as.character(iteration), color = model_run))+
ggplot2::scale_color_manual(values = c("#D65F00", "black", "blue")) +
ggplot2::scale_linetype_manual(values = rep("solid", 50)) +
ggplot2::guides(linetype = FALSE) +
ggplot2::facet_wrap(. ~ scenario) +
Expand Down
94 changes: 17 additions & 77 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,70 +149,18 @@ dir.create(run_SSMSE_dir)

The cod model with H = 0.65 (as in scenario 1) is included as external
package data in SSMSE. However, we will need to modify it to use as an
operating model with H = 1 (as in scenario 2).
operating model with H = 1 (as in scenario 2). Note in this case that
refit\_OM is false, so the model is not being refit, just run through
without fitting. To condition the new model on the same data as the
input model, refit\_OM should be TRUE.

``` r
cod_mod_path <- system.file("extdata", "models", "cod", package = "SSMSE")
# copy to a new location:
file.copy(from = cod_mod_path, to = run_SSMSE_dir, recursive = TRUE)
## [1] TRUE
file.rename(from = file.path(run_SSMSE_dir, "cod"), to = file.path(run_SSMSE_dir, "cod-1"))
## [1] TRUE
cod_1_path <- file.path(run_SSMSE_dir, "cod-1")
# make model read initial values from control file and not ss.par
start <- r4ss::SS_readstarter(file = file.path(cod_1_path, "starter.ss"), verbose = FALSE)
start$init_values_src # verify reading from the control file
## [1] 0
# change the steepness from 0.2 to 0.1 in the control files
r4ss::SS_changepars(dir = cod_1_path, ctlfile = "control.ss_new",
newctlfile = "control_modified.ss", strings = "SR_BH_steep", newvals = 1)
## parameter names in control file matching input vector 'strings' (n=1):
## [1] "SR_BH_steep"
## These are the ctl file lines as they currently exist:
## LO HI INIT PRIOR PR_SD PR_type PHASE env_var&link dev_link dev_minyr
## 108 0.2 1 0.65 0.7 0.05 0 -4 0 0 0
## dev_maxyr dev_PH Block Block_Fxn Label Linenum
## 108 0 0 0 0 SR_BH_steep 108
## line numbers in control file (n=1):
## 108
## wrote new file to control_modified.ss with the following changes:
## oldvals newvals oldphase newphase oldlos newlos oldhis newhis oldprior
## 1 0.65 1 -4 -4 0.2 0.2 1 1 0.7
## newprior oldprsd newprsd oldprtype newprtype comment
## 1 0.7 0.05 0.05 0 0 # SR_BH_steep
# read in the par file to save the historic recruitment deviations:
parfile <- r4ss::SS_readpar_3.30(parfile = file.path(cod_1_path, "ss.par"),
datsource = file.path(cod_1_path, "ss3.dat"),
ctlsource = file.path(cod_1_path, "control.ss_new"),
verbose = FALSE)
# remove files with old steepness values
file.remove(file.path(cod_1_path, "control.ss_new"))
## [1] TRUE
file.remove(file.path(cod_1_path, "control.ss"))
## [1] TRUE
file.remove(file.path(cod_1_path, "ss.par"))
## [1] TRUE
file.rename(from = file.path(cod_1_path, "control_modified.ss"),
to = file.path(cod_1_path, "control.ss"))
## [1] TRUE
```

Rerun this model with no estimation to get valid ss.par and
control.ss\_new files, then add in the historical recruitment
deviations:

``` r
SSMSE:::run_ss_model(dir = cod_1_path,
admb_options = "-maxfn 0 -phase 50 -nohess",
verbose = FALSE)
# add back original recdevs into the model (b/c not specified through the ctl file)
new_parfile <- r4ss::SS_readpar_3.30(parfile = file.path(cod_1_path, "ss.par"),
datsource = file.path(cod_1_path, "ss3.dat"),
ctlsource = file.path(cod_1_path, "control.ss"), verbose = FALSE)
# add in the recdevs to new the parfile
new_parfile$recdev1[, "recdev"] <- parfile$recdev1[, "recdev"]
r4ss::SS_writepar_3.30(new_parfile, outfile = file.path(cod_1_path, "ss.par"),
verbose = FALSE)
develop_OMs(OM_name = "cod", out_dir = run_SSMSE_dir, par_name = "SR_BH_steep",
par_vals = 1, refit_OMs = FALSE, hess = FALSE)
# OM model for scenario 2
cod_1_path <- file.path(run_SSMSE_dir, "cod_SR_BH_steep_1")
```

### Examine the management procedure used
Expand Down Expand Up @@ -382,7 +330,7 @@ run_SSMSE(scen_name_vec = c("H-ctl", "H-1"),# name of the scenario
scope = "2", # to use the same recruitment devs across scenarios.
impl_error_pattern = "none", # Don't use implementation error
run_EM_last_yr = FALSE, # Run the EM in 106
run_parallel = TRUE,
run_parallel = FALSE,
sample_struct_list = sample_struct_list, # How to sample data for running the EM.
seed = 12345) #Set a fixed integer seed that allows replication
```
Expand Down Expand Up @@ -436,31 +384,23 @@ library(dplyr)
```

``` r
summary$ts <- tidyr::separate(summary$ts,
col = model_run,
into = c(NA, "model_type"),
remove = FALSE,
sep = "_",
extra = "drop")
# check values for cod_OM
summary$scalar %>%
dplyr::filter(iteration == 1) %>%
dplyr::filter(scenario == "H-1") %>%
dplyr::select(iteration, scenario, SR_BH_steep, model_run)
## iteration scenario SR_BH_steep model_run
## 1 1 H-1 1.00 cod-1_OM
## 2 1 H-1 0.65 cod_EM_103
## 3 1 H-1 0.65 cod_EM_106
## 4 1 H-1 0.65 cod_EM_init
## iteration scenario SR_BH_steep model_run
## 1 1 H-1 0.65 cod_EM_103
## 2 1 H-1 0.65 cod_EM_init
## 3 1 H-1 1.00 cod_SR_BH_steep_1_OM


# plot SSB by year and model run - need to correct using code from the
# think tank
ggplot2::ggplot(data = subset(summary$ts, model_run %in% c("cod_OM", "cod-1_OM", "cod_EM_106")),
# plot SSB by year and model run
ggplot2::ggplot(data = subset(summary$ts, model_run %in% c("cod_OM", "cod_SR_BH_steep_1_OM", "cod_EM_103")),
ggplot2::aes(x = year, y = SpawnBio)) +
ggplot2::geom_vline(xintercept = 100, color = "gray") +
ggplot2::geom_line(ggplot2::aes(linetype = as.character(iteration), color = model_type))+
ggplot2::scale_color_manual(values = c("#D65F00", "black")) +
ggplot2::geom_line(ggplot2::aes(linetype = as.character(iteration), color = model_run))+
ggplot2::scale_color_manual(values = c("#D65F00", "black", "blue")) +
ggplot2::scale_linetype_manual(values = rep("solid", 50)) +
ggplot2::guides(linetype = FALSE) +
ggplot2::facet_wrap(. ~ scenario) +
Expand Down
Binary file modified man/figures/README-plot_SSB-1.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified man/figures/README-plot_rel_SSB-1.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

0 comments on commit ee336f5

Please sign in to comment.