#' Fit the \code{aspic} biomass dynamic model.
#'
#' @description
#' Fits the aspic model to catch and catch per unit effort data
#'
#' @param object; an \code{aspic} object
#' @param dir; an optional \code{dir} where aspic text files used for fitting can be found
#' @return An aspic object with fitted values and parameter estimates
#' @seealso \code{\link{aspic},\link{biodyn},\link{boot},\link{jk}}
#' @export
#' @examples
#' ## get aspic input file (i.e. *.inp) for North Atlantic swordfish
#' file=paste(system.file(package="aspic"),"exdata/swon.inp",sep="/")
#'
#' swon=aspic(file)
#'
#' swon=fit(swon)
#'
#' plot(swon)
setMethod('fit', signature(object='aspic',index="missing"),
function(object, dir=tempdir(), package=class(object), exeNm="aspic",jk=FALSE,copyExe=FALSE)
runExe(object=object, dir=dir, package=package, exeNm=exeNm,jk=jk,copyExe=copyExe))
#' jk, jack knifes \code{aspic}
#'
#' @description
#' Fits the aspic model to catch and catch per unit effort data removing 1 cpue observation at a time
#'
#' @param object; an \code{aspic} object or
#' @param object; a character string giving an aspic "inp" file
#' @param dir; an optional \code{dir} where aspic text files used for fitting can be found
#' @return An aspic object with fitted values and parameter estimates
#' @seealso \code{\link{biodyn},\link{boot},\link{fit}}
#' @export
#' @rdname profile
#' @examples
#' \dontrun{
#' data(asp)
#' asp=jk(asp)}
setMethod('jk', signature(object='aspic'),
function(object, dir=tempdir(), package=class(object), exeNm="aspic")
runExe(object=object, dir=dir, package=package, exeNm=exeNm,jk=TRUE))
#' boot, Bootstraps the ASPIC biomass dynamic model.
#'
#' @description
#' Bootstraps the aspic model
#'
#' @param object; an \code{aspic} object or
#' @param object; a character string giving an aspic "inp" file
#' @param dir; an optional \code{dir} where aspic text files used for fitting can be found
#' @return An aspic object with fitted values and parameter estimates
#' @seealso \code{\link{biodyn},\link{boot},\link{jk}}
#' @export
#' @rdname boot
#' @examples
#' \dontrun{
#' data(asp)
#' asp=boot(asp)}
setMethod('boot', signature(object='aspic'),
function(object, dir=tempdir(), package=class(object), exeNm="aspic",boot=500)
runBoot(object=object, dir=dir,package=package, exeNm=exeNm,boot=boot))
utils::globalVariables(c("year","swon","year","B","obs"))
chkIters=function(object){
N=max(dims(object)$iter,dims(object@control)$iter,dims(object@params)$iter)
#params(object)=FLPar(object@control[,"val",drop=T])
if (N>1){
object@stock=propagate(stock( object),N)
if (dims(object@params)[1]!=N)
object@params=propagate(object@params,N)
object@ll =propagate(object@ll, N)
object@objFn=propagate(object@objFn,N)
}
object@stock=propagate(stock(object),dims(object)$iter)
object@stock=window(stock(object),end=max(as.numeric(dimnames(catch(object))$year))+1)
return(object)}
jkIdx=function(x) dimnames(x)[[1]][ !is.na(x$index)]
runExe=function(object,package="aspic",exeNm=package,dir=tempdir(),jk=FALSE,copyExe=FALSE){
object@index=object@index[object@index$year %in% range(object)["minyear"]:range(object)["maxyear"],]
if (any(is.na(object@catch))){
tmp=ddply(object@index, .(year), with, sum(catch,na.rm=TRUE))
object@catch=as.FLQuant(tmp[,"V1"], dimnames=list(year=tmp[,"year"]))
dmns=dimnames(object@catch)
dmns$year=c(dmns$year,as.numeric(max(dmns$year))+1)
object@stock=FLQuant(NA,dimnames=dmns)
}
#oldwd =setExe(exeNm,package,dir)
#oldwd =setExe(exeNm,package,dir)
oldwd=getwd()
setwd(dir)
path=biodyn:::exe("aspic")
if (.Platform$OS.type == "windows" & copyExe)
file.copy(paste(paste(system.file("bin", "windows", package=package, mustWork=TRUE),exeNm, sep="/"),"exe",sep="."), dir)
else if (copyExe)
file.copy( paste(system.file("bin", "linux", package=package, mustWork=TRUE),exeNm, sep="/"), dir)
## Jack knife if wished
j=1
if (jk){
object=propagate(object,length(jkIdx(object@index)))
object@params=propagate(object@params,length(jkIdx(object@index)))
j = jkIdx(object@index)
index=object@index}
us=paste("u",seq(length(dimnames(params(object))$params[grep("q",dimnames(params(object))$params)])),sep="")
dmns=list(params=us,quantity=c("ll","ss","n"),iter=seq(1))
object@ll=FLPar(array(NA,dim=unlist(lapply(dmns,length)),dimnames=dmns))
object=chkIters(object)
for (i in seq(dims(object)$iter)){
m_ply(c("prn","rdat","bio","inp","fit","sum","rdatb","det","sum","bot"), function(x)
if (file.exists(paste(exeNm,".",x,sep=""))) system(paste("rm ",exeNm,".",x,sep="")))
if (jk){
object@index=index
object@index[j[i],"index"]=NA
}
# create exe input files
.writeAspicInp(FLCore:::iter(object,i),what="FIT",niter=1,fl=paste(exeNm,".inp",sep=""))
# run
#system(paste("./", exeNm, paste(" ",exeNm,".inp",sep=""),sep=""))
#print(paste("aspic", paste(" ","aspic",".inp",sep=""),sep=""))
system(paste("aspic", paste(" ","aspic",".inp",sep=""),sep=""))
#system(paste(exeNm, paste(" ",exeNm,".inp",sep=""),sep=""))
rdat=dget(paste(exeNm,"rdat",sep="."))
#rdat$estimates
object@params[c("b0","msy","k"),i]=rdat$estimates[c("B1.K","MSY","K")]
object@params[4:dim(object@params)[1],i]=rdat$estimates[substr(names(rdat$estimates),1,2)=="q."]
names(rdat$t.series)=tolower(names(rdat$t.series))
FLCore:::iter(object@stock,i)=as.FLQuant(transform(rdat$t.series[,c("year","b")],data=b)[c("year","data")])[,dimnames(object@stock)$year]
if (.Platform$OS!="windows"){
try(object@objFn[1,i]<-rdat$diagnostics$obj.fn.value)
#try(object@objFn[2,i]<-rdat$diagnostics$rsquare)
rtn=try(aspicPrn(paste(exeNm,"prn",sep=".")))
if (is.data.frame(rtn)) object@diags=rtn[!is.na(rtn$residual),]
object@diags=transform(object@diags,stock. =hat/c(object@params[grep("q",dimnames(params(object))$params),i])[name],
stockHat=obs/c(object@params[grep("q",dimnames(params(object))$params),i])[name])
object@diags=merge(object@diags,model.frame(mcf(FLQuants(stock =object@stock,
harvest=harvest(object))),drop=TRUE),all=T)
object@diags$name=ac(factor(object@diags$name, labels=unique(object@index$name)))
object@diags$stock=object@diags$stock.
object@diags=object@diags[,-10]
object@diags=object@diags[!is.na(object@diags$name),]
} else {
rtn=try(readAspic(paste(exeNm,"prn",sep=".")))
if (is.data.frame(rtn)) object@diags=rtn[!is.na(rtn$residual),]
object@diags=transform(object@diags,stock. =hat/c(object@params[grep("q",dimnames(params(object))$params),i])[name],
stockHat=obs/c(object@params[grep("q",dimnames(params(object))$params),i])[name])
object@diags$name=ac(factor(object@diags$name, labels=unique(object@index$name)))
object@diags=merge(object@diags,model.frame(mcf(FLQuants(stock=object@stock,harvest=harvest(object))),drop=TRUE),all=T)
#object@diags$stock=object@diags$stock.
object@diags=object@diags[,-10]
try(object@objFn[1,i]<-sum(diags(object)$residual^2,na.rm=T))
object@diags=object@diags[!is.na(object@diags$name),]
}
dgs=subset(object@diags,!is.na(object@diags$residual))
try(object@ll@.Data[,"ll",i]<-daply(dgs, .(name), with, biodyn:::calcLogLik(residual,residual/residual,type=3)))
try(object@ll@.Data[,"ss",i]<-daply(dgs, .(name), with, sum(residual^2)))
try(object@ll@.Data[,"n", i]<-daply(dgs, .(name), function(x) dim(x)[1]))
}
#if (dims(object)$iter!=1) object@diags=data.frame(NULL)
setwd(oldwd)
return(object)}
runBoot=function(object, package="aspic", exeNm=package, dir=tempdir(),boot=500){
object@index=object@index[object@index$year %in% range(object)["minyear"]:range(object)["maxyear"],]
if (boot<3) {
boot=3
warning("Requires a minimum of 3 bootstraps so boot option changed")
}
## add catch baed on index catches if missing
if (any(is.na(object@catch))){
tmp=ddply(object@index, .(year), with, mean(catch,na.rm=TRUE))
object@catch=as.FLQuant(tmp[,"V1"], dimnames=list(year=tmp[,"year"]))
}
## add catch baed on index catches if missing
if (dim(object@control)[3] >1) stop("control can only have iter dim of 1")
if (dim(object@params)[2]==1) object@params=propagate(object@params,boot)
if (dim(object@params)[2]!=boot) stop("params iters either have to be 1 or same as number of boot")
if (dims(object@catch)$iter>1) stop("catch must only have 1 iter")
dmns=dimnames(object@catch)
dmns$year=c(dmns$year,as.numeric(max(dmns$year))+1)
object@stock=propagate(FLQuant(NA,dimnames=dmns),boot)
oldwd=biodyn:::setExe(package,exeNm,dir)
m_ply(c("prn","rdat","bio","inp","fit","sum","rdatb","det","sum","bot"), function(x)
if (file.exists(paste(exeNm,".",x,sep=""))) system(paste("rm ",exeNm,".",x,sep="")))
# create exe input files
.writeAspicInp(object,what="BOT",niter=boot,fl=paste(exeNm,".inp",sep=""))
# run
system(paste("./", exeNm, paste(" ",exeNm,".inp",sep=""),sep=""))
det=aspicDet(paste(exeNm,"det",sep="."))
bio=aspicBio(paste(exeNm,"bio",sep="."))
coerceDP=function(x) FLPar(unlist(c(t(x))),params=names(x),iter=dim(x)[1])
parNms=c("b0",biodyn:::modelParams(tolower(model(object))))
object@params[parNms,]=coerceDP(det[,parNms])
qNms=names(det)[!(names(det) %in% c("iter","stock","harvest","r","trial","loss","msy","bmsy","brel","frel","b1.k",parNms))]
object@params[-(seq(length(parNms)))][]=unlist(c(det[,qNms]))
object@stock=bio[["stock"]]%*%bio[["bmsy"]]
setwd(oldwd)
return(object)}
setMethod('fit', signature(object='aspics',index="missing"),
function(object, dir=tempdir(), package=class(object), exeNm="aspic",jk=FALSE,
.combine=NULL,
.multicombine=T,.maxcombine=10,.packages=c("aspic","plyr")){
if (is.null(.combine)) ..combine=list else ..combine=.combine
res=foreach(i=names(object), .combine=..combine,
.multicombine=.multicombine,
.maxcombine =.maxcombine,
.packages =.packages) %dopar% {
wkdir=tempfile('file', dir)
dir.create(wkdir, showWarnings = FALSE)
fit(object[[i]],dir=wkdir)}
if (is.null(.combine)) {
res=aspics(res)
names(res)=names(object)}
res})
setMethod('boot', signature(object='aspics'),
function(object, dir=tempdir(), package=class(object), exeNm="aspic",boot=500,
.combine=NULL,
.multicombine=T,.maxcombine=10,.packages=c("aspic","plyr")){
if (is.null(.combine)) ..combine=list else ..combine=.combine
res=foreach(i=names(object), .combine=..combine,
.multicombine=.multicombine,
.maxcombine =.maxcombine,
.packages =.packages) %dopar% {
wkdir=tempfile('file', dir)
dir.create(wkdir, showWarnings = FALSE)
boot(object[[i]],dir=wkdir,boot=boot)}
if (is.null(.combine)) {
res=aspics(res)
names(res)=names(object)}
res})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.