utils::globalVariables(c("aspicCtl","k","iUAspic","value","."))
utils::globalVariables('read.table')
utils::globalVariables('read.table')
utils::globalVariables('qqnorm')
checkFile=function(x){
ln=tolower(scan(x,nlines=1,what=character(),sep="\n"))
if (any(maply(data.frame(pattern=c("trial","loss","msy","bmsy","brel","frel"),stringsAsFactors=F), grep, x=ln)>0))
return("det")
}
#setMethod('writeAspic', signature(object="aspic"),
# function(object,index=object@index,what="FIT",niter=1,fl="aspic.inp",...) .writeAspicInp(object,index,what,niter,fl=fl,...))
.writeAspicInp<-function(object,index=object@index,what="FIT",niter=ifelse(what=="FIT",1,501),fl="aspic.inp"){
dgts=options()$digits
options(digits=22)
if (file.exists("aspic.fit")) file.remove("aspic.fit")
if (file.exists("aspic.prn")) file.remove("aspic.prn")
if (file.exists("aspic.rdat")) file.remove("aspic.rdat")
if (file.exists("aspic.sum")) file.remove("aspic.sum")
dmmy=expand.grid(year=min(as.numeric(as.character(index$year))):max(as.numeric(as.character(index$year))),
name=unique(index$name))[,2:1]
u=merge(dmmy,index,all=TRUE,sort=FALSE)
u$index[is.na(u$index)]=-9999
u$catch[is.na(u$catch)]=0
comment=rep("",22)
comment[ 1]= "\n"
comment[ 2]= "\n"
comment[ 3]= "\n"
comment[ 4]= "\t512 ## Verbosity\n"
comment[ 5]= "\t## Number of bootstrap trials, <= 1000\n"
comment[ 6]= "\t## 0=no MC search, 1=search, 2=repeated srch; N trials\n"
comment[ 7]= "\t## Convergence crit. for simplex\n"
comment[ 8]= "\t## Convergence crit. for restarts, N restarts\n"
comment[ 9]= "\t## Conv. crit. for F; N steps/yr for gen. model\n"
comment[10]= "\t## Maximum F when cond. on yield\n"
comment[11]= "\t## Stat weight for B1>K as residual (usually 0 or 1)\n"
comment[12]= "\t## Number of fisheries (data series)\n"
comment[13]= "\t## Statistical weights for data series\n"
comment[14]= "\t## B1/K (starting guess, usually 0 to 1)\n"
comment[15]= "\t## MSY (starting guess)\n"
comment[16]= "\t## K (carrying capacity) (starting guess)\n"
comment[17]= "\t## q (starting guesses -- 1 per data series)\n"
comment[18]= "\t## Estimate flags (0 or 1) (B1/K,MSY,K,q1...qn)\n"
comment[19]= "\t## Min and max constraints -- MSY\n"
comment[20]= "\t## Min and max constraints -- K\n"
comment[21]= "\t## Random number seed\n"
comment[22]= "\t## Number of years of data in each series\n"
# search trials simplex restarts nrestarts effort nsteps maxf
#1e+00 1e+05 1e-08 3e-08 6e+00 1e-04 0e+00 8e+00
cat(what ,comment[ 1],file=fl,append=FALSE)
cat("FLR generated" ,comment[ 2],file=fl,append=TRUE)
cat(ac(object@model), ac(object@conditioning), ac(object@obj) ,comment[ 3],file=fl,append=TRUE)
cat( comment[ 4],file=fl,append=TRUE)
cat(niter ,comment[ 5],file=fl,append=TRUE)
cat(as.integer(object@options[c("search","trials")]) ,comment[ 6],file=fl,append=TRUE)
cat(object@options["simplex"] ,comment[ 7],file=fl,append=TRUE)
cat(object@options["restarts"],as.integer(object@options["nrestarts"]),comment[ 8],file=fl,append=TRUE)
cat(object@options["effort"],object@options["nsteps"],comment[ 9],file=fl,append=TRUE)
cat(object@options["maxf"] ,comment[10],file=fl,append=TRUE)
cat(0 ,comment[11],file=fl,append=TRUE)
cat(dim(object@control)[1]-3 ,comment[12],file=fl,append=TRUE)
cat(object@control[-(1:3),"lambda"] ,comment[13],file=fl,append=TRUE)
cat(object@control["b0", "val"] ,comment[14],file=fl,append=TRUE)
cat(object@control["msy", "val"] ,comment[15],file=fl,append=TRUE)
cat(object@control["k", "val"] ,comment[16],file=fl,append=TRUE)
cat(object@control[-(1:3),"val"] ,comment[17],file=fl,append=TRUE)
cat(object@control[ ,"fit"] ,comment[18],file=fl,append=TRUE)
cat(object@control["msy",c("min","max")] ,comment[19],file=fl,append=TRUE)
cat(object@control["k", c("min","max")] ,comment[20],file=fl,append=TRUE)
cat(as.integer(object@rnd) ,comment[21],file=fl,append=TRUE)
cat(daply(u,.(name), with, length(name)),comment[22],file=fl,append=TRUE)
d_ply(u,.(name), function(x) {
if(any(names(x)=="type"))
names(x)[seq(length(names(x)))[names(x)=="type"]]="code"
cat(as.character(unique(x$name)),"\n",file=fl,append=TRUE)
cat(as.character(x$code[1]),"\n",sep="",file=fl,append=TRUE)
cat(apply(x[,c("year","index","catch")],1,paste, collapse=" "),sep="\n",file=fl,append=TRUE)})
# l_ply(index,
# function(idx){
# cat(name(idx) ,"\n",file=fl,append=TRUE)
# cat(type(idx) ,"\n",file=fl,append=TRUE)
#
# mm=model.frame(FLQuants(col2=index(idx),col3=catch.n(idx)),drop=T)
# mm[is.na(mm)]=-1
#
# cat(paste(t(apply(as.matrix(mm),1,paste,collapse=" ")),"\n"),file=fl,append=TRUE)})
options(digits=dgts)
return()}
## local function to calculated expected QQ line
qqLine <- function(x,y){
qtlx <- quantile(x, probs=c(0.25,0.75), na.rm=T)
qtly <- quantile(y, probs=c(0.25,0.75), na.rm=T)
a <- (qtly[1]- qtly[2]) / (qtlx[1] - qtlx[2])
b <- qtly[1] - qtlx[1] * a
res <- c(a,b)
names(res) <- NULL
names(res) <- c("a","b")
return(res)}
fnDiags=function(res){
res$residualLag <- c(res$residual[-1],NA)
qq. <- qqnorm(res$residual,plot.it=FALSE,na.rm=T)
res$qqx <- qq.$x
res$qqy <- qq.$y
qqpar <- qqLine(qq.$x,qq.$y)[c("a","b")]
res$qqHat=qqpar["a"]*res$qqx+qqpar["b"]
res}
files=data.frame(ext =c("inp","bio","prb","rdat","det","prn","fit"),
desc=c("Input file with data, starting guesses, and run settings",
"Estimated B and F trajectory for each bootstrap trial",
"As .bio but with projection results",
"Inputs and estimates specially formatted for R",
"Estimates from each bootstrap trial",
"index fitted and observed",
"results in print format"),stringsAsFactor=FALSE)
getExt <- function(file)
tolower(substr(file,max(gregexpr("\\.", file)[[1]])+1,nchar(file)))
checkExt=function(x) (tolower(getExt(x)) %in% files[,"ext"])
#### Aspic #####################################################################################
.readAspic<-function(x,relative=TRUE){
if (!checkExt(x)) stop(cat("File", x, "does not have a valid extension"))
type=getExt(x)
return(switch(tolower(type),
"inp" =aspicInp(x),
"bio" =aspicBio(x,relative),
"prb" =aspicPrb(x,relative),
"rdat"=aspicRdat(x),
"ctl" =aspicCtl(x),
"det" =aspicDet(x),
"prn" =aspicPrn(x),
"fit" =aspicFit(x)))
}
################################################################################
aspicBio =function(file,relative=TRUE){
t. <-scan(file,skip=4)
nits<-scan(file,skip=1,nmax=1)
yrs <-scan(file,skip=2,nmax=2)
nyrs<-diff(yrs)
nval<-nyrs*2+3
yrs <-yrs[1]:yrs[2]
b. <-FLQuant(t.[unlist(tapply(((1:nits)-1)*nval+2, 1:nits,function(x,y=nyrs+1) x:(x+y-1)))],dimnames=list(year=yrs, iter=1:nits))
f. <-FLQuant(t.[unlist(tapply(((1:nits)-1)*nval+nyrs+4,1:nits,function(x,y=nyrs) x:(x+y-1)))],dimnames=list(year=yrs[-length(yrs)], iter=1:nits))
bmsy<-FLQuant(t.[unlist(tapply(((1:nits)-1)*nval+1, 1:nits,function(x,y=1) x:(x+y-1)))],dimnames=list( iter=1:nits))
fmsy<-FLQuant(t.[unlist(tapply(((1:nits)-1)*nval+nyrs+3,1:nits,function(x,y=1) x:(x+y-1)))],dimnames=list( iter=1:nits))
if (relative){
b.=b.%/%bmsy
f.=f.%/%fmsy}
return(FLQuants(stock=b.,harvest=f.,bmsy=bmsy,fmsy=fmsy))}
aspicPrb =function(file,relative=TRUE){
## Stuff
nits<-scan(file,skip=1,nmax=1)
yrs <-scan(file,skip=2,nmax=2)
t. <-scan(file,skip=4)
ncol<-yrs[2]-yrs[1]+2
## stock
first<-rep((1:nits-1)*ncol*2,each=yrs[2]-yrs[1]+1)+(1:(ncol-1))+1
b. <-FLQuant(t.[first],dimnames=list(year=yrs[1]:yrs[2],iter=1:nits))
if (relative){
first<-((1:nits-1)*ncol*2)+1
bmsy <-FLQuant(t.[first],dimnames=list(iter=1:nits))
b. <-sweep(b.,6,bmsy,"/")
}
## F
first<-rep((1:nits-1)*ncol*2+ncol,each=yrs[2]-yrs[1]+1)+(1:(ncol-1))+1
f. <-FLQuant(t.[first],dimnames=list(year=yrs[1]:yrs[2],iter=1:nits))[,ac(yrs[1]:(yrs[2]-1))]
if (relative)
{
first<-((1:nits-1)*ncol*2)+ncol+1
fmsy <-FLQuant(t.[first],dimnames=list(iter=1:nits))
f. <-sweep(f.,6,fmsy,"/")
}
return(FLQuants(harvest=FLQuant(f.),
stock =FLQuant(b.)))}
aspicRdat=function(file){
res=dget(file)
names(res$estimates)=tolower(names(res$estimates))
names(res$estimates)[seq(length(res$estimates))[names(res$estimates)=="b.bmsy"]]="bbmsy"
names(res$estimates)[seq(length(res$estimates))[names(res$estimates)=="f.fmsy"]]="ffmsy"
names(res$estimates)[seq(length(res$estimates))[names(res$estimates)=="yield.eq"]]="yieldeq"
names(res$estimates)[seq(length(res$estimates))[names(res$estimates)=="b1.k"]] ="b0"
names(res$t.series)=tolower(names(res$t.series))
names(res$t.series)[seq(length(res$t.series))[names(res$t.series)=="f.total"]] ="harvest"
names(res$t.series)[seq(length(res$t.series))[names(res$t.series)=="b"]] ="stock"
names(res$t.series)[seq(length(res$t.series))[names(res$t.series)=="l.tot.obs"]]="catch"
names(res$t.series)[seq(length(res$t.series))[names(res$t.series)=="f.fmsy"]] ="ffmsy"
names(res$t.series)[seq(length(res$t.series))[names(res$t.series)=="b.bmsy"]] ="bbmsy"
return(res)}
aspicDet =function(x){
det=read.table(x,header=TRUE)
nms=names(det)
names(det)[seq(length(nms))[nms=="trial"]]="iter"
names(det)[seq(length(nms))[nms=="b1.k"]]="b0"
names(det)[seq(length(nms))[nms=="brel"]]="stock"
names(det)[seq(length(nms))[nms=="frel"]]="harvest"
if (all(c("msy","k") %in% nms) & !("r" %in% nms)) det=transform(det,r=4*msy/k)
dim(det)[2]
det=det[,c(1:8,dim(det)[2],10:dim(det)[2]-1)]
dNms=names(det)
dNms=gsub("q\\.0","q",dNms)
dNms=gsub("q\\.","q",dNms)
names(det)=dNms
det}
aspicInp =function(x){
res=aspic()
aspicC=function(file) {
inp=scan(file,sep="\n",what=character())
inp=str_trim(sub("\t+"," ",inp))
ctrl=mlply(inp[5:21], function(x) {
tmp=strsplit(x," ")
tmp=unlist(tmp)[nchar(unlist(tmp))>0]})
ctrl=llply(ctrl,as.numeric)
ctrl=llply(ctrl,function(x) x[!is.na(x)])
return(ctrl)}
ctrl=suppressWarnings(aspicC(x))
ops =strsplit(scan(x,sep="\n",what=character(),skip=2,nlines=1), " ")[[1]]
ops =ops[nchar(ops)>0]
res@model[1] =factor(ops[1])
res@conditioning=factor(ops[2])
res@obj =factor(ops[3])
# [8] "7 ## Number of fisheries (data series)"
n =ctrl[[8]]
params=FLPar("b0"=as.numeric(NA),"k"=as.numeric(NA),"msy"=as.numeric(NA))
parNms=c(c("b0","msy","k"),paste("q",seq(n),sep=""))
res@params=FLPar(as.numeric(NA),parNms,iter=1)
res@control=FLPar(array(as.numeric(NA),c(length(c(c("b0","msy","k"),paste("q",seq(n),sep=""))),5,1),dimnames=list(params=parNms,c("fit","min","val","max","lambda"),iter=1)))
# [10] "1.00000 ## B1/K (starting guess, usually 0 to 1)"
res@control["b0", "val"]=ctrl[[10]][1]
# [11] "3.0000E+04 ## MSY (starting guess)"
res@control["msy","val"]=ctrl[[11]]
# [12] "2.6700E+05 ## K (carrying capacity) (starting guess)"
res@control["k", "val"]=ctrl[[12]]
# [13] "2.1126E-06 6.0195E-06 9.7627E-06 1.4944E-04 2.9980E-06 4.2138E-04 8.3406E-04 ## q (starting guesses -- 1 per data series)"
res@control[parNms[-(1:3)],"val"]=ctrl[[13]][1:n]
# [14] "0 1 1 1 1 1 1 1 1 1 ## Estimate flags (0 or 1) (B1/K,MSY,K,q1...qn)"
res@control[,"fit"]=ctrl[[14]]
# [15] "1.0000E+02 1.0000E+07 ## Min and max constraints -- MSY"
res@control["msy",c("min","max")]=ctrl[[15]]
# [16] "1.0000E+04 2.0000E+07 ## Min and max constraints -- K"
res@control["k", c("min","max")]=ctrl[[16]]
res@control[parNms[-(1:3)],"min"]=res@control[parNms[-(1:3)],"val"]*0.01
res@control[parNms[-(1:3)],"max"]=res@control[parNms[-(1:3)],"val"]*100
res@control["b0","min"]=0.01
res@control["b0","max"]=1
# [9] "1.0000E+00 1.0000E+00 1.0000E+00 1.0000E+00 1.0000E-02 1.0000E+00 1.0000E-02 ## Statistical weights for data series"
if (length(ctrl[[9]])==0)
res@control[parNms[-(1:3)],"lambda"][]=1.0 else
res@control[parNms[-(1:3)],"lambda"]=ctrl[[9]]
# [17] "6745260 ## Random number seed
res@rnd =ctrl[[17]]
#
# # [3] "1.0000E-08 ## Convergence crit. for simplex"
# if (length(ctrl[[3]])==0) ctrl[[3]]=1.0-08
# res@conv["simplex"]=ctrl[[3]]
#
# # [4] "3.0000E-08 6 ## Convergence crit. for restarts, N restarts"
# res@conv["restart"]=ctrl[[4]][1]
# res@nRestart=ctrl[[4]][1]
#
# # [5] "1.0000E-04 0 ## Conv. crit. for F; N steps/yr for gen. model"
# res@conv["F"]=ctrl[[5]][1]
# res@nSteps =ctrl[[5]][1]
#
# # [6] "8.0000 ## Maximum F when cond. on yield"
# res@maxF =ctrl[[6]]
# # [7] "0.0 ## Stat weight for B1>K as residual (usually 0 or 1)"
# res@wt =ctrl[[7]][1]
res@index=iUAspic(x)#,"aspic")
res@catch=as.FLQuant(ddply(index(res),.(year), with,
data.frame(data=sum(catch))))
res@stock= FLQuant(NA,dimnames=dimnames(res@catch))
res@stock=window(res@stock,end=dims(res@catch)$maxyear+1)
rng =range(res@index$year)
names(rng)=c("minyear","maxyear")
res@range =rng
return(res)}
aspicPrn =function(x){
#x="/home/lkell/Dropbox/MyStuff/WHM/analysis/Inputs/aspic/Base case runs/whmrun1bb.prn"
res=read.table(x,header=TRUE,colClasses="numeric")
res=res[,seq(dim(res)[2]-2)]
obs=melt(res[,seq((dim(res)[2]-1)/2+1)],id.var="year")
est=melt(res[,c(1,((dim(res)[2]-1)/2+2):dim(res)[2])],id.var="year")
res=data.frame(transform(obs,obs=value,name=gsub(".obs","",obs$variable))[,c("year","name","obs")],
hat=est$value)
if (is.factor(res$hat)) res$hat=as.numeric(as.character(res$hat))
res$residual=log(res$obs/res$hat)
res=ddply(res,.(name),fnDiags)
names(res)[3]="obs"
res}
#x="/home/laurie/Desktop/flr/tests/aspic/Inputs/swon/2009/run9/aspic.fit"
aspicFit=function(x){
txt=str_trim(scan(x,sep="\n",what=as.character()))
start=seq(length(txt))[substr(txt,1,9)=="ESTIMATED"]+5
end =seq(length(txt))[substr(txt,1,7)=="RESULTS"] -3
txt =txt[start:end]
res=as.data.frame(t(matrix(as.numeric(unlist(strsplit(txt," +"))),ncol=end-start+1,nrow=10)))[,-1]
names(res)=c("year","harvest","biomass","biomassMn","yield","yieldHat","sp","harvestMSY","biomassMSY")
res}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.