mpb: Management Procedures

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}