Skip to content

Commit 583a5c2

Browse files
Ovec8hkinOvec8hkin
committed
New plotting code in ggplot
1 parent 9dc4b45 commit 583a5c2

17 files changed

Lines changed: 1307 additions & 1074 deletions

‎functions/fun_read_dat.R‎

Lines changed: 287 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,287 @@
1+
# fun_read_dat.r
2+
# Created by John Trochta
3+
# Date created: 07/27/2020
4+
# Summary:
5+
# Project population dynamics to next year with given catches
6+
7+
library(tidyverse)
8+
9+
# Read BASA data files (.dat and .ctl) into list format for easy parameter access.
10+
read.data.files <- function(dat.dir){
11+
12+
filename <- vector(length=5)
13+
filename[1] <- paste0(dat.dir, "PWS_ASA.dat")
14+
filename[2] <- paste0(dat.dir, "PWS_ASA(ESS).ctl")
15+
filename[3] <- paste0(dat.dir, "PWS_ASA(covariate).ctl")
16+
filename[4] <- paste0(dat.dir, "agecomp_samp_sizes.txt")
17+
filename[5] <- paste0(dat.dir, "PWS_ASA_disease.dat")
18+
19+
PWS_ASA.dat <- data_reader(filename=filename[1])
20+
PWS_ASA_ESS.ctl <- data_reader(filename=filename[2])
21+
PWS_ASA_covariate.ctl <- data_reader(filename=filename[3])
22+
agecomp_samp_sizes.txt <- data_reader(filename=filename[4])
23+
PWS_ASA_disease.dat <- data_reader(filename=filename[5])
24+
25+
par.names <- c("nyr", "nyr_tobefit", "nage", "waa", "fecundity",
26+
"pound_catch", "pk", "foodbait_catch", "gillnet_catch", "seine_yield", "perc.female",
27+
"mdm", "egg", "egg_se", "adfg_hydro_year_start", "adfg_hydro", "pwssc_hydro_year_start", "pwssc_hydro", "pwssc_hydro_se",
28+
"seine_age_comp", "spawn_age_comp", "juvenile_survey")
29+
names(PWS_ASA.dat) <- par.names
30+
31+
par.names <- c("ctl_file", "seine_ess", "spawn_ess", "vhsv_ess", "ich_ess")
32+
names(PWS_ASA_ESS.ctl) <- par.names
33+
34+
par.names <- c("std_covs", "num_recruit_cov", "recruit_fixed_rand", "cov_on", "regime_shift_89", "r_beta_change",
35+
"num_mort_covs", "mort_fixed_rand", "mort_on", "mort_age_impact", "disease_covs", "winter_mort_devs", "m_btea_change")
36+
names(PWS_ASA_covariate.ctl) <- par.names
37+
38+
par.names <- c("seine_sample_size", "spawn_sample_size", "vhsv_sample_size", "ich_smaple_size")
39+
names(agecomp_samp_sizes.txt) <- par.names
40+
41+
par.names <- c("vhsv_age_prevalence", "vhsv_obs_start", "vhsv_est_start", "vhsv_recov_prob",
42+
"ich_age_prevalence", "ich_obs_start", "ich_est_start", "ich_recov_prob")
43+
names(PWS_ASA_disease.dat) <- par.names
44+
45+
return(listN(PWS_ASA.dat,
46+
PWS_ASA_ESS.ctl,
47+
PWS_ASA_covariate.ctl,
48+
agecomp_samp_sizes.txt,
49+
PWS_ASA_disease.dat))
50+
}
51+
52+
# Read and coerce BASA data into appropriate R data type (vector, list, matrix, etc.).
53+
data_reader <- function(filename) {
54+
# The user needs to make sure there is a blank line at the end of
55+
# each file and that each data type (vector number or matrix) is
56+
# separated by a blank line
57+
58+
# This is kind of convoluted
59+
text <- readLines(filename)
60+
values <- grep("^\\s{0,2}[0-9]",text)
61+
signed.values <- grep("^\\s{0,2}[-]",text)
62+
read.these <- sort(c(values,signed.values))
63+
nlines <- length(text)
64+
indices <- seq(1:nlines)
65+
indices <- indices[read.these]
66+
first.differences <- c(diff(indices),5) # This accounts for the last data
67+
data.types <- length(first.differences[first.differences>1])
68+
69+
data <- vector("list", data.types)
70+
j <- 1
71+
temp <- NA
72+
for(i in 1:length(indices)){
73+
temp.1 <- scan(filename,skip=indices[i]-1,nlines=1,quiet=TRUE,flush=FALSE)
74+
if(first.differences[i]>1){
75+
if(all(is.na(temp))){
76+
data[[j]] <- temp.1
77+
j <- j+1
78+
} else{
79+
temp <- rbind(temp,temp.1)
80+
data[[j]] <- temp
81+
rownames(data[[j]]) <- NULL
82+
temp <- NA
83+
j <- j+1
84+
}
85+
} else{
86+
if(all(is.na(temp))){
87+
temp <- temp.1
88+
} else{
89+
temp <- rbind(temp,temp.1)
90+
}
91+
}
92+
}
93+
return(data)
94+
}
95+
96+
# Read BASA parameter estimates (from .par file) into name list
97+
read.par.file <- function(filename){
98+
library(stringr)
99+
par.vals <- data_reader(filename)
100+
text <- readLines(filename)
101+
par.names <- text[grep("^#", text)]
102+
par.names <- apply(as.array(par.names[2:length(par.names)]), 1, function(x) stringr::str_match(x, "[a-zA-Z0-9_]+"))
103+
names(par.vals) <- par.names
104+
return(par.vals)
105+
}
106+
107+
read.biomass.estimates <- function(model.dir, nyr=NA){
108+
fname <- paste0(here::here(model.dir), "/mcmc_out/PFRBiomass.csv")
109+
biomass.data <- read.table(fname, header = FALSE, sep = ",", dec=".")
110+
111+
nyr <- ifelse(is.na(nyr), ncol(biomass.data)-1, nyr)
112+
years <- 1980:(1980+nyr)
113+
114+
colnames(biomass.data) <- years
115+
116+
return(biomass.data[,1:nyr])
117+
118+
}
119+
120+
read.exploit.rates <- function(model.dir, nyr=NA){
121+
raw.data <- read.data.files(model.dir)
122+
123+
pfrb <- read.biomass.estimates(model.dir, nyr=nyr)
124+
nyr <- ifelse(is.na(nyr), ncol(pfrb), nyr)
125+
126+
total.catch.biomass <- compute.catch.biomass(raw.data$PWS_ASA.dat, nyr)
127+
exploit.rate <- t(total.catch.biomass/t(pfrb))
128+
129+
return(exploit.rate)
130+
131+
}
132+
133+
compute.catch.biomass <- function(data, nyr){
134+
135+
weight.at.age <- data$waa[1:nyr,]
136+
137+
fb.nya <- data$foodbait_catch[1:nyr,]
138+
pound.nya <- data$pound_catch[1:nyr,]
139+
gillnet.nya <- data$gillnet_catch[1:nyr,]
140+
seine.yield <- data$seine_yield[1:nyr]
141+
142+
fb.nya.biomass <- weight.at.age * fb.nya
143+
pound.nya.biomass <- weight.at.age * pound.nya
144+
gillnet.nya.biomass <- weight.at.age * gillnet.nya
145+
146+
fb.biomass.annual <- rowSums(fb.nya.biomass) # Now sum the biomass over all age classes for each year
147+
pound.biomass.annual <- rowSums(pound.nya.biomass)
148+
gillnet.biomass.annual <- rowSums(gillnet.nya.biomass)
149+
150+
fb.biomass.annual <- replace(fb.biomass.annual, fb.biomass.annual == 0, NA)
151+
pound.biomass.annual <- replace(pound.biomass.annual, pound.biomass.annual == 0, NA)
152+
gillnet.biomass.annual <- replace(gillnet.biomass.annual, gillnet.biomass.annual == 0, NA)
153+
seine.yield <- replace(seine.yield, seine.yield == 0, NA)
154+
155+
# Matrix of catches by gear type in mt
156+
total.catch <- cbind(fb.biomass.annual, pound.biomass.annual, gillnet.biomass.annual, seine.yield)
157+
total.catch[is.na(total.catch)] <- 0
158+
total.catch.biomass <- rowSums(total.catch) # total catches by year in mt
159+
160+
return(total.catch.biomass)
161+
162+
}
163+
164+
read.survey.estimates <- function(model.dir){
165+
166+
nage=10
167+
168+
fnames <- list(
169+
mdm = paste0(model.dir, "/MDM.csv"),
170+
pwssc.hydro = paste0(model.dir, "/HYD_PWSSC.csv"),
171+
adfg.hydro = paste0(model.dir, "/HYD_ADFG.csv"),
172+
spawn.age.comp = paste0(model.dir, "/SpAC.csv"),
173+
seine.age.comp = paste0(model.dir, "/SeAC.csv"),
174+
egg = paste0(model.dir, "/EGG.csv"),
175+
juv.schools = paste0(model.dir, "/juv_schools.csv")
176+
)
177+
178+
spac.data <- read.csv(fnames$spawn.age.comp, header=FALSE)
179+
spac.data <- apply(spac.data, 2, median)
180+
spac.data <- matrix(spac.data, ncol=nage, byrow=TRUE)
181+
182+
seac.data <- read.csv(fnames$seine.age.comp, header=FALSE)
183+
seac.data <- apply(seac.data, 2, median)
184+
seac.data <- matrix(seac.data, ncol=nage, byrow=TRUE)
185+
186+
mdm.data <- read.csv(fnames$mdm, header=FALSE)
187+
mdm.data <- apply(mdm.data, 2, median)
188+
mdm.data <- as.vector(mdm.data)
189+
190+
pwssc.hydro.data <- read.csv(fnames$pwssc.hydro, header=FALSE)
191+
pwssc.hydro.data <- apply(pwssc.hydro.data, 2, median)
192+
pwssc.hydro.data <- as.vector(pwssc.hydro.data)
193+
194+
adfg.hydro.data <- read.csv(fnames$adfg.hydro, header=FALSE)
195+
adfg.hydro.data <- apply(adfg.hydro.data, 2, median)
196+
adfg.hydro.data <- as.vector(adfg.hydro.data)
197+
198+
egg.data <- read.csv(fnames$egg, header=FALSE)
199+
egg.data <- apply(egg.data, 2, median)
200+
egg.data <- as.vector(egg.data)
201+
egg.data[egg.data == 0] <- -9
202+
203+
juv.schools.data <- read.csv(fnames$juv.schools, header=FALSE)
204+
juv.schools.data <- apply(juv.schools.data, 2, median)
205+
juv.schools.data <- as.vector(juv.schools.data)
206+
207+
return(listN(mdm.data, egg.data, pwssc.hydro.data, adfg.hydro.data, juv.schools.data, spac.data, seac.data))
208+
}
209+
210+
read.catch.data <- function(cr, sims, nyr){
211+
212+
fnames <- c("foodbait_catch.csv", "gillnet_catch.csv", "pound_catch.csv", "seine_catch.csv")
213+
214+
data <- data.frame(year=NA, catch=NA, fishery=NA, control.rule=NA, sim=NA)
215+
for(s in sims){
216+
for(f in fnames){
217+
fname <- paste0(here::here("results/"), cr, "/sim_", s, "/year_", nyr, "/results/", f)
218+
dat.fname <- paste0(here::here("results/"), cr, "/sim_", s, "/year_", nyr, "/model/")
219+
220+
waa <- read.data.files(dat.fname)$PWS_ASA.dat[[4]]
221+
waa <- waa[(42+1):(42+1+nyr-1),]
222+
catch.data <- as.matrix(read.csv(fname))[1:nyr, 2:11]
223+
catch.biomass <- apply(catch.data*waa, 1, sum)
224+
225+
d <- data.frame(year=1:nyr, catch=catch.biomass, fishery=str_split(f, "_")[[1]][1], control.rule=cr, sim=s)
226+
data <- data %>% bind_rows(d)
227+
}
228+
}
229+
return(data)
230+
231+
}
232+
233+
accumulate.results.data <- function(nyr.sim, total.sims, seeds, trial, fnames, byage=FALSE, ncols=1){
234+
235+
if(length(byage) == 1){
236+
byage <- rep(byage, length(fnames))
237+
}
238+
239+
#d.vec <- ifelse(byage, 10, 1)
240+
241+
data.matrices <- vector("list", length(fnames))
242+
for(i in 1:length(fnames)){
243+
d <- ifelse(byage[i], 10, ncols)
244+
data.matrices[[i]] <- array(NA, dim=c(nyr.sim, d, total.sims))
245+
}
246+
247+
for(i in 1:length(fnames)){
248+
for(s in 1:length(seeds)){
249+
f <- paste0(here::here("results"), "/", trial, "/sim_", seeds[s], "/year_", nyr.sim, "/results/", fnames[i])
250+
data.matrices[[i]][,,s] <- as.matrix(read.csv(f)[1:nyr.sim,-1])
251+
}
252+
253+
}
254+
255+
return(data.matrices)
256+
257+
}
258+
259+
accumulate.assessment.posteriors <- function(nyr.sim, total.sims, seeds, trial, fname="PFRBiomass.csv"){
260+
261+
f1 <- paste0(here::here("results"), "/", trial, "/sim_", seeds[1], "/year_1/model/mcmc_out/", fname)
262+
n.samps <- nrow(read.csv(f1, header=FALSE))
263+
#print(n.samps)
264+
data.matrices <- vector("list", length(fname))
265+
for(i in 1:length(fname)){
266+
data.matrices[[i]] <- array(NA, dim=c(n.samps, nyr.sim, total.sims))
267+
}
268+
269+
for(i in 1:nyr.sim){
270+
for(s in 1:length(seeds)){
271+
f <- paste0(here::here("results"), "/", trial, "/sim_", seeds[s], "/year_", i-1, "/model/mcmc_out/", fname)
272+
print(f)
273+
data <- read.csv(f, header=FALSE)
274+
data.matrices[[1]][,i,s] <- as.matrix(data[(nrow(data)-n.samps+1):nrow(data), ncol(data)])
275+
}
276+
}
277+
278+
return(data.matrices)
279+
}
280+
281+
282+
# Function for simultaneously creating names of variables/elements within a list
283+
listN <- function(...){
284+
anonList <- list(...)
285+
names(anonList) <- as.character(substitute(list(...)))[-1]
286+
anonList
287+
}

‎plotting/compute_plot_products.R‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
library(here)
2+
library(tidyverse)
3+
source(paste0(here::here("functions/"), "fun_read_dat.R"))
4+
5+
6+
7+
compute.exploit.rate <- function(model.dir, nyr, years){
8+
exploit.rate <- read.exploit.rates(model.dir, nyr)
9+
10+
exploit.rate.df <- as_tibble(exploit.rate) %>%
11+
pivot_longer(everything(), names_to="year", values_to="exploit") %>%
12+
group_by(year) %>%
13+
median_qi(exploit, .width=c(0.95)) %>%
14+
print(n=10)
15+
16+
exploit.zeros <- exploit.rate.df[exploit.rate.df$exploit == 0, ]
17+
18+
return(listN(exploit.rate.df, exploit.zeros))
19+
}
20+
21+
compute.biomass.traj <- function(model.dir, nyr, years){
22+
biomass <- read.biomass.estimates(model.dir, nyr)
23+
24+
biomass.df <- as_tibble(biomass) %>%
25+
pivot_longer(everything(), names_to="year", values_to="biomass") %>%
26+
group_by(year) %>%
27+
median_qi(biomass, .width=c(0.10, 0.50, 0.95)) %>%
28+
print(n=10)
29+
30+
prob.below.threshhold <- apply(biomass, 2, function(x) sum(x < 19958))/nrow(biomass)
31+
32+
biomass.df$prob <- as.vector(rep(prob.below.threshhold, 3))
33+
34+
return(biomass.df)
35+
}
36+
37+
compute.pfrb.posterior <- function(model.dir, nyr, years){
38+
biomass <- read.table(paste0(model.dir, "mcmc_out/PFRBiomass.csv"), header = FALSE, sep = ",", dec=".")[,nyr]
39+
40+
biomass.df <- data.frame(biomass=biomass)
41+
prob.below.threshold <- round(sum(biomass < 20000)/length(biomass), 2)
42+
43+
biomass.quants <- round(as.vector(apply(as.matrix(biomass), 2, quantile, c(0.025, 0.5, 0.975)))/1000, 2)
44+
45+
return(listN(biomass.df, biomass.quants, prob.below.threshold))
46+
}
47+
48+
compute.recruitment <- function(model.dir, nyr, years){
49+
age.3.recruits <- read.table(paste0(model.dir, "mcmc_out/Age3.csv"), header = FALSE, sep = ",", dec=".")[,1:nyr]
50+
colnames(age.3.recruits) <- years
51+
52+
age.3.recruits.df <- as_tibble(age.3.recruits) %>%
53+
pivot_longer(everything(), names_to="year", values_to="recruits") %>%
54+
group_by(year) %>%
55+
median_qi(recruits, .width=c(0.50, 0.95)) %>%
56+
print(n=10)
57+
58+
return(age.3.recruits.df)
59+
}
60+
61+
compute.catch.biomass <- function(model.dir, nyr, years){
62+
63+
data <- read.data.files(model.dir)
64+
65+
weight.at.age <- data$waa[1:nyr,]
66+
67+
fb.nya <- data$foodbait_catch[1:nyr,]
68+
pound.nya <- data$pound_catch[1:nyr,]
69+
gillnet.nya <- data$gillnet_catch[1:nyr,]
70+
seine.yield <- data$seine_yield[1:nyr]
71+
72+
fb.nya.biomass <- weight.at.age * fb.nya
73+
pound.nya.biomass <- weight.at.age * pound.nya
74+
gillnet.nya.biomass <- weight.at.age * gillnet.nya
75+
76+
fb.biomass.annual <- rowSums(fb.nya.biomass) # Now sum the biomass over all age classes for each year
77+
pound.biomass.annual <- rowSums(pound.nya.biomass)
78+
gillnet.biomass.annual <- rowSums(gillnet.nya.biomass)
79+
80+
fb.biomass.annual <- replace(fb.biomass.annual, fb.biomass.annual == 0, NA)
81+
pound.biomass.annual <- replace(pound.biomass.annual, pound.biomass.annual == 0, NA)
82+
gillnet.biomass.annual <- replace(gillnet.biomass.annual, gillnet.biomass.annual == 0, NA)
83+
seine.yield <- replace(seine.yield, seine.yield == 0, NA)
84+
85+
# Matrix of catches by gear type in mt
86+
total.catch <- cbind(fb.biomass.annual, pound.biomass.annual, gillnet.biomass.annual, seine.yield)
87+
total.catch[is.na(total.catch)] <- 0
88+
total.catch.biomass <- rowSums(total.catch) # total catches by year in mt
89+
90+
return(total.catch.biomass)
91+
92+
}

0 commit comments

Comments
 (0)