#' @title read_result_hake
#'
#' @description Reads the ADMB report (result) file from the Namibian hake assessment
#'
#' @export
#'
#' @param infilename A character string incuding the path (directory of the input file name.
#' For hake the name of the file is normally namhakdata.dat
#' @param outfilename A character string, including path (directory) to the report
#' file. The report file has normally the ending .rep.
#' @param runname A character specifying the name of the run
read_rbx_hake <- function(infilename,outfilename,runname) {
txt <- readLines(outfilename)
txt <- stringr::str_trim(txt)
tmpfile <- tempfile()
# SigmaR_input SigmaR_output MSY Bsp/K
p1 <- grep("SigmaR_input",txt)
x <- str_clean(txt,p1,p1+1,header=T)
x <- reshape2::melt(x)
names(x) <- c("V1","V2")
# Alpha_first
p1 <- grep("Alpha_first",txt)
p2 <- p1
g <- txt[(p1):(p2)]
g <- stringr::str_trim(g)
g <- unlist(strsplit(g," "))
x2_names <- g[c(1,3,5,7,9)]
x2_values <- as.numeric(g[c(2,4,6,8,10)])
x2 <- data.frame(V1=x2_names,V2=x2_values)
sr_parameters <- x2
names(sr_parameters) <- c("variable","value")
x <- rbind(x,x2)
# 'ln(Likelihoods)
p1 <- grep("Likelihoods",txt) + 1
p2 <- grep("M ",txt)
x2 <- str_clean(txt,p1,p2)
likelihoods <- x2[1:9,]
names(likelihoods) <- c("variable","value")
refpts <- x2[c(10,11,15,16,17),]
names(refpts) <- c("variable","value")
refpts <- rbind(refpts,data.frame(variable="Fmsy",value=refpts$value[5]/refpts$value[4]))
refpts$details <- c("Virgin ssb","Virgin exploitable biomass",
"SSB at MSY","Exploitable biomass at MSY",
"Maximum sustainable yield",
"Fishing mortality resulting in MSY")
x <- rbind(x,x2)
p1 <- grep("sigmaCPUE1",txt)
p2 <- grep("Addvariance",txt)
x2 <- str_clean(txt,p1,p2)
x <- rbind(x,x2)
p1 <- grep("sigCAA_com",txt)
p2 <- grep("sigSeal_Index",txt)
x2 <- str_clean(txt,p1,p2)
x <- rbind(x,x2)
p1 <- grep("Natural mortality by age",txt) + 1
p2 <- grep("M\\(8\\)",txt)
x2 <- str_clean(txt,p1,p2,header=F)
x <- rbind(x,x2)
p1 <- grep("CPUEseries1",txt)
p2 <- grep("Survseries2",txt)
x2 <- str_clean(txt,p1,p2,header=F)
x2$V1 <- paste(x2$V1,x2$V2,sep="_")
x2$V4 <- paste(x2$V1,x2$V4,sep="_")
x2a <- x2[,c(1,3)]
names(x2a) <- c("V1","V2")
x <- rbind(x,x2a)
x2a <- x2[,c(4,5)]
names(x2a) <- c("V1","V2")
x <- rbind(x,x2a)
x <- x[!duplicated(x),]
names(x) <- c("variable","value")
tmp_out <- x
# INDICES - the first results by year
p1 <- grep('Year obs1 pred1 obs2 pred2 obs3 pred3 obs4 pred4 obs5',txt)
p2 <- grep('Virgin Age Structure',txt)
g <- txt[(p1):(p2-2)]
# amend the first line (header)
g[1] <- paste('year','oU1','pU1','oU2','pU2','oU3','pU3','oU4','pU4','oU5','pU5',
'oU6','pU6','oUs','pUs','oUw','pUw','oUss','xx','pUss',sep=' ')
rby <- str_clean(g,1,length(g),sep=' ',header=T)
# The second last colum (now labelled 'xx' is redundant). Lets get rid of it:
rby$xx <- NULL
# ZERO's are actually NA - get rid of them
rby <- reshape2::melt(rby,id.vars='year')
rby$value[rby$value == 0] <- NA
rby <- reshape2::dcast(rby,year ~ variable,value.var='value')
# STOCK IN NUMBERS
p1 <- grep('Numbers at Age',txt) + 1
p2 <- grep('Year Bsp Bexp F RecRes Catch Depletion TotalB',txt) - 2
rbya <- str_clean(txt,p1,p2,header=F)
names(rbya) <- c('year',0:8)
rbya <- reshape2::melt(rbya,id.vars='year')
names(rbya) <- c('year','age','n')
rbya$age <- as.integer(as.character(rbya$age))
# SOME MORE results by year
p1 <- grep('Year Bsp Bexp F RecRes Catch Depletion TotalB',txt)
p2 <- grep('Commercial selectivity',txt)
tmp <- str_clean(txt,p1,p2-2,header=T)
names(tmp) <- c('year','ssb','biofish','f','rR','y','dep','bio')
rby <- plyr::join(rby,tmp,by="year")
# COMMERCIAL SELECTION PATTERN
p1 <- p2 + 1
p2 <- grep("Observed catch in numbers by age",txt) - 2
sel <- str_clean(txt,p1,p2,header=T,fill=T)
names(sel) <- c('year',0:8)
sel <- reshape2::melt(sel,id.vars='year')
names(sel) <- c('year','age','sel')
sel$age <- as.numeric(as.character(sel$age))
rbya <- plyr::join(rbya,sel,by=c("year","age"))
# HERE WILL SKIP READING SOME FILES - BECAUSE SAME INFO COMES LATER IN THE OUTPUT
p1 <- grep('Replacement Yield by year',txt)
p2 <- grep('----------Stock-recruitment curve',txt)
x <- str_clean(txt,p1+1,p2-1,header=F)
names(x) <- c('year','rY','bio','landings') # some guesswork
x <- x[,c("year","rY","landings")]
rby <- plyr::join(rby,x,by="year")
# Add recruitment to rby
rby$rec <- rbya[rbya$age == 0,"n"]
# Reorder column - just for aesthetic reasons
rby <- rby[,c("year","rec","bio","ssb","biofish","f","landings","rY","dep",
"oU1","pU1","oU2","pU2","oU3","pU3","oU4","pU4","oU5","pU5",
"oU6","pU6","oUs","pUs","oUw","pUw","oUss","pUss")]
rby$landings[rby$landings == 0] <- NA
rby$f[rby$f == 0] <- NA
# Create a separate object for the observed vs fitted value
x <- rby[,c(1,10:ncol(rby))]
x.obs <- x[,c(1,2,4,6,8,10,12,14,16,18)]
x.pre <- x[,c(1,3,5,7,9,11,13,15,17,19)]
cn <- c("year","ICSEAF (1.3 + 1.4)","ICSEAF -1.5","GLM","winter Spanish",
"summer Spanish","7vessel CPUE","summer","winter","seal scat")
names(x.obs) <- names(x.pre) <- cn
x.obs <- reshape2::melt(x.obs,"year")
names(x.obs) <- c("year","fleet","obs")
x.pre <- reshape2::melt(x.pre,"year")
names(x.pre) <- c("year","fleet","pre")
rby.fit <- plyr::join(x.obs,x.pre,by=c("year","fleet"))
# STOCK RECRUITMENT STUFF
p1 <- grep('----------Stock-recruitment curve',txt)
p2 <- grep('----------CommercialCAA',txt)
xy_ssb.r <- str_clean(txt,p1+1,p2-2,header=F)
xy_ssb.r <- as.data.frame(t(xy_ssb.r))
names(xy_ssb.r) <- c('ssb','r')
#### Catch at age
p1 <- grep('----------CommercialCAA',txt)
p2 <- grep('----------Survey1',txt)[2]
x <- str_clean(txt,p1+1,p2-1)
names(x) <- c('year','age','oC','pC','rC')
rbya <- plyr::join(rbya,x,by=c("year","age"))
# Summer survey at age
p1 <- grep('----------Survey1',txt)[2]
p2 <- grep('----------Survey2',txt)[2]
x <- str_clean(txt,p1+1,p2-1,header=F)
names(x) <- c('year','age','oU1','pU1','rU1')
rbya <- plyr::join(rbya,x,by=c("year","age"))
#### Winter survey residuals
p1 <- grep('----------Survey2',txt)[2]
p2 <- grep('biomass of fish by age',txt)
x <- str_clean(txt,p1+1,p2-2,header=F)
names(x) <- c('year','age','oU2','pU2','rU2')
rbya <- plyr::join(rbya,x,by=c("year","age"))
#### Biomass at age
p1 <- grep('biomass of fish by age',txt)
p2 <- grep('biomass of fish that died due to natural mortality by age for J-P Roux',txt)
x <- str_clean(txt,p1+1,p2-2,header=F)
g <- txt[(p1+1):(p2-2)]
writeLines(g,tmpfile)
x <- read.table(tmpfile,header=FALSE)
names(x) <- c('year',0:8)
x <- reshape2::melt(x,id.vars='year')
names(x) <- c('year','age','b')
x$age <- as.integer(as.character(x$age))
rbya <- plyr::join(rbya,x,by=c("year","age"))
if(!missing(runname)) {
rby$name <- runname
rbya$name <- runname
xy_ssb.r$name <- runname
tmp_out$Name <- runname
}
## Extract some results by age data
# Natural mortality
i <- tmp_out$variable %in% c("M(0)","M(1)","M(2)","M(3)","M(4)","M(5)","M(6)","M(7)","M(8)")
rba <- tmp_out[i,]
rba$age <- rep(0:8)
rba <- rba[,c("age","value")]
names(rba) <- c("age","m")
rbya$m <- rep(rba$m,(max(rbya$year)-min(rbya$year)+1))
# quick fix here
x <- read_input_hake(infilename)
names(x) <- c("year","age","sW","cW","mat")
rbya <- plyr::join(rbya,x,by=c("year","age"))
i <- x$year == min(x$year)
x <- x[i,]
rba$sW <- x$sW
rba$cW <- x$cW
rba$mat <- x$mat
rbya$yc <- rbya$year - rbya$age
rownames(rba) <- NULL
rbx <- list(rby=rby,rba=rba,rby.fit=rby.fit,rbya=rbya,xy_ssb.r=xy_ssb.r,
likelihoods=likelihoods,refpts=refpts,tmp_out=tmp_out)
return(rbx)
}
#' @title read_input_hake
#'
#' @description Reads the ADMB input data file from the Namibian hake assessment
#'
#' @export
#'
#' @param filename A character string, including path (directory) to the report
#' file. The report file has normally the ending .rep.
read_input_hake <- function(filename) {
tmpfile <- tempfile()
txt <- readLines(filename)
txt <- stringr::str_trim(txt)
# Weight-at-age
# Start-Year
# 0 1 2 3 4 5 6 7 8
p1 <- grep('Weight-at-age',txt)
p2 <- grep('Mid-Year',txt)
g <- txt[(p1+3):(p2-1)]
g <- stringr::str_trim(g)
g <- stringr::str_replace_all(g,"\\s+","\t")
#g <- as.data.frame(g)
writeLines(g,tmpfile)
w1 <- as.data.frame(read.table(tmpfile,header=FALSE))
txt <- readLines(filename)
txt <- stringr::str_trim(txt)
# Weight-at-age
# Start-Year
# 0 1 2 3 4 5 6 7 8
p1 <- grep('Mid-Year',txt)
p2 <- grep('Catches',txt)
g <- txt[(p1+2):(p2-1)]
g <- stringr::str_trim(g)
g <- stringr::str_replace_all(g,"\\s+","\t")
#g <- as.data.frame(g)
writeLines(g,tmpfile)
w2 <- as.data.frame(read.table(tmpfile,header=FALSE))
names(w1) <- names(w2) <- c("year",0:8)
d <- reshape2::melt(w1,c("year"))
names(d) <- c("year","age","w1")
w2 <- reshape2::melt(w2,c("year"))
d$w2 <- w2$value
d$age <- as.integer(as.character(d$age))
# Maturity
p1 <- grep('#\tEnd\tof\tfile\tflag',txt) - 2
p2 <- p1 + 1
x <- str_clean(txt,p1,p2,header=T)
x <- reshape2::melt(x)
names(x) <- c("age","mat")
x$age <- as.integer(as.character(x$age))
d <- plyr::join(d,x,by="age")
return(d)
}
#' @title read_sim_hake
#'
#' @description Reads the simulation files for Hake.
#'
#' @return Returns a list containing two data.frames, the raw data (raw) and
#' summarized data (qua). The latter contains the 0.05, 0.25, 0.50, 0.75 and
#' 0.95 percentiles as well as the mean.
#'
#' @export
#'
#' @param path The directory path where the Files are located files are located
#'
read_sim_hake <- function(path) {
# fingerplots
filename <- dir(path,pattern="fingerplots",full.names = TRUE)
p1 <- stringr::str_locate(filename,"fingerplots_")[,2] + 1
p2 <- stringr::str_locate(filename,".out")[,1] - 1
runname <- as.numeric(stringr::str_sub(filename,p1,p2))
for (i in 1:length(filename)) {
x <- read.table(filename[i])
names(x) <- c("year","TAC","B_B1990","exp")
n <- length(unique(x$year))
iter <- nrow(x)/n
x$iter <- rep(1:iter,each=n)
x$run <- runname[i]
if(i == 1) {
raw <- x
} else {
raw <- rbind(raw,x)}
}
# Economic
filename <- dir(path,pattern="Economic",full.names = TRUE)
p1 <- stringr::str_locate(filename,"Economic_")[,2] + 1
p2 <- stringr::str_locate(filename,".out")[,1] - 1
runname <- as.numeric(stringr::str_sub(filename,p1,p2))
for (i in 1:length(filename)) {
x <- read.table(filename[i])
names(x) <- c("year","TAC","value","employm","profit","vessel1","vessel2")
n <- length(unique(x$year))
iter <- nrow(x)/n
x$iter <- rep(1:iter,each=n)
x$run <- runname[i]
if(i == 1) {
raw2 <- x
} else {
raw2 <- rbind(raw2,x)}
}
raw <- plyr::join(raw,raw2[,-2],by=c("year","iter","run"))
## Depletion
filename <- dir(path,pattern="Depletion",full.names = TRUE)
p1 <- stringr::str_locate(filename,"Depletion_")[,2] + 1
p2 <- stringr::str_locate(filename,".out")[,1] - 1
runname <- as.numeric(stringr::str_sub(filename,p1,p2))
for (i in 1:length(filename)) {
x <- read.table(filename[i])
names(x) <- rep(c("year","TAC","depletion"),ncol(x)/3)
j <- seq(1,ncol(x),by=3)
year <- tac <- dep <- vector()
for (k in 1:length(j)) {
year <- c(year,x[,j[k]])
tac <- c(tac,x[,j[k]+1])
dep <- c(dep,x[,j[k]+2])
}
x <- data.frame(year=year,depletion=dep)
n <- length(unique(x$year))
iter <- nrow(x)/n
x$iter <- rep(1:iter,n)
x$run <- runname[i]
if(i == 1) {
raw2 <- x
} else {
raw2 <- rbind(raw2,x)}
}
raw <- plyr::join(raw,raw2,by=c("year","iter","run"))
q <- reshape2::melt(raw,c("year","iter","run"))
q <- plyr::ddply(q,c("year","run","variable"),summarize,
q05=quantile(value,0.05,na.rm=T),
q25=quantile(value,0.25,na.rm=T),
q50=quantile(value,0.50,na.rm=T),
q75=quantile(value,0.75,na.rm=T),
q95=quantile(value,0.95,na.rm=T),
m=mean(value,na.rm=T))
return(list(raw=raw,qua=q))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.