-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrunEstGrowth.R
179 lines (151 loc) · 6.47 KB
/
runEstGrowth.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
###############################################################################
###############################################################################
#-----------------------------------------------------------------------------#
####Author : Kelli Faye Johnson
####Contact : [email protected]
####Lastupdate :
####Purpose :
####Packages :
####Inputs :
####Outputs :
####Remarks : Character width = 80
#-----------------------------------------------------------------------------#
###############################################################################
###############################################################################
###############################################################################
## Step 01
## Variable inputs: objects in Step 01 may need alteration prior to running
###############################################################################
#For Christine's machine
main.path <- "c:/users/christine stawitz/documents/github"
dir.main <- file.path(main.path,"estgrowth")
dir.dropbox <- "c:/users/christine stawitz/documents/dropbox/ExtGrowthResults"
#For Kelli's machine
## Set the working directory
if (Sys.info()["user"] == "kelli") {
dir.main <- "c:/ss/estgrowth"
dir.dropbox <- "c:/users/kelli/dropbox/estgrowth"
}
if (Sys.info()["user"] == "kfjohns") {
dir.main <- "t:/estgrowth"
dir.dropbox <- "t:/estgrowth"
}
# Which species do you want to use?
my.spp <- c("hake", "flatfish", "yellow")
my.spp <- my.spp[order(my.spp)]
my.totnum <- 1:150
# Logical whether or not to run ss3sim bias adjustment to gain information
# about recruitment from years with more information
my.bias.num <- 10
# Register the number of cores you want to use
my.corenum <- ifelse(Sys.info()["user"] == "kfjohns", 10, 3)
# Logical whether in testing mode or not
testingmode <- FALSE
###############################################################################
## Step 02
## Load packages
###############################################################################
library(devtools)
install_github("r4ss/r4ss", ref = "master")
install_github("ss3sim/ss3sim", ref = "master")
install_github("ss3sim/ss3models", ref = "master")
library(doParallel); library(foreach); library(r4ss);
library(ss3models); library(ss3sim);
dir.models <- system.file("models", package = "ss3models")
dir.sub <- file.path(dir.main, "test")
dir.cases <- file.path(dir.main, "cases")
dir.results <- file.path(dir.main, "results")
# Copy case files from github growth_models repository to the local folder
# for casefiles.
setwd(dir.main)
file.copy(system.file("cases", package = "ss3models"), ".", recursive = TRUE)
###############################################################################
## Step 03
## Model locations
###############################################################################
# file below need true values of M for each species
truem <- vector(length = length(my.spp))
for (ind in seq_along(my.spp)) {
om <- file.path(ss3model(my.spp[ind], "om"), "ss3.ctl")
pars <- SS_parlines(om)
truem[ind] <- pars[grep("NatM", pars$Label), "INIT"]
}
source(file.path(dir.main, "lib", "createcasefiles.r"))
my.casefiles <- list(A = "agecomp", C = "calcomp", D = "mlacomp",
E = "E", F = "F",
I = "index", L = "lcomp") #,
# Below are cases that we are not currently using
# S == selectivity case files (dome)
# M == time-varying natural mortality
# S = c(toupper(rev(letters)[1:6])), M = "M")
# Import scenario list from excel file
setwd("lib")
done <- system("make")
if (done == 127) stop("User must manually update lib\\scenario.csv")
setwd("..")
scenariosfile <- file.path(dir.main, "lib", "scenarios.csv")
scenarios <- read.table(scenariosfile, sep = ",", header = TRUE)
torun <- unlist(lapply(scenarios$complete, paste0, my.spp))
#set working directory
dir.create(dir.sub, showWarnings = FALSE)
setwd(dir.sub)
if (testingmode) {
# # # Run a single iteration of a test scenario
test <- paste("A31-C0-D0-E1-F0-I0-L31", my.spp, sep = "-")
unlink(test, recursive = TRUE)
for (ind in seq_along(my.spp)) {
use.om <- ss3model(my.spp[ind], "om")
use.em <- ss3model(my.spp[ind], "em")
run_ss3sim(iterations = 50, scenarios = test[ind],
case_folder = dir.cases, case_files = my.casefiles,
om_dir = use.om, em_dir = use.em, bias_adjust = FALSE,
ignore.stdout = TRUE, show.output.on.console = FALSE)
}
}
if (length(ls(pattern = "test$", envir=.GlobalEnv)) > 0) {
removeme <- file.path(dir.sub, test)
ignore <- sapply(removeme, unlink, recursive = TRUE)
}
# Set up running in parallel if specified
registerDoParallel(cores = my.corenum)
getDoParWorkers() # check
# Determine if some iterations have already been run
# If so only run more iterations (to the maximum my.totrun)
# of the F1 fishing scenario.
done <- length(dir(dir(pattern = my.spp[1], full.names = TRUE)[1],
pattern = "[0-9]$"))
if(done > 0) {
my.totnum <- seq(done + 1, max(my.totnum), by = 1)
torun <- torun[grepl("F1", torun)]
}
#Use the following to run all combinations
for(s in seq_along(my.spp)){
#finds the appropriate folder for each species
use.om <- ss3model(my.spp[s], "om")
use.em <- ss3model(my.spp[s], "em")
#run scenarios that include the given species
use.scen <- torun[sapply(my.spp[s], grepl, torun)]
run_ss3sim(iterations = my.totnum, scenarios = use.scen,
case_folder = dir.cases, case_files = my.casefiles,
om_dir = use.om, em_dir = use.em,
bias_already_run = ifelse(any(grepl("^1$", my.totnum)), FALSE, TRUE),
bias_adjust = ifelse(my.bias.num > 0 & any(grepl("^1$", my.totnum)), TRUE, FALSE),
bias_nsim = my.bias.num,
ignore.stdout = TRUE, show.output.on.console = FALSE,
parallel = ifelse(getDoParWorkers() > 1, TRUE, FALSE))
# Should also maybe set ss_mode = "optimized"
}
###############################################################################
## Step
## Get results
###############################################################################
get_results_all(overwrite_files = FALSE,
parallel = ifelse(getDoParWorkers() > 1, TRUE, FALSE))
scalars <- read.csv(dir(pattern = "r.csv", full.names = TRUE))
ts <- read.csv(dir(pattern = "s.csv", full.names = TRUE))
file.copy(dir(pattern = "r.csv", full.names = TRUE),
file.path(dir.dropbox, "scalars.csv"), overwrite = TRUE, copy.mode = TRUE)
file.copy(dir(pattern = "scalar", full.names = TRUE),
file.path(dir.dropbox, "ts.csv"), overwrite = TRUE, copy.mode = TRUE)
source(file.path(dir.main, "lib", "results.r"))
setwd(dir.main)