utils::globalVariables(c('ldply','melt','variable'))
globalVariables("lambda")
#if (!isGeneric('control'))
setGeneric('control', function(object,method,...) standardGeneric('control'))
#if (!isGeneric('control<-'))
setGeneric('control<-', function(object,value) standardGeneric('control<-'))
#' setParams<-
#'
#' Add parameters for catchability and sigma in the param slot given an index or indices
#' of stock abundance as \code{FLQuant} or \code{FLQuants}.
#'
#' @param object an object of class \code{biodyn}
#' @param value \code{params}
#'
#' @export
#' @rdname setParams
#'
#' @aliases setParams<-,biodyn,FLPar-method setParams<-,biodyn,FLQuant-method setParams<-,biodyn,FLQuants-method setParams<-,biodyn,data.frame-method setParams<- setParams<-,biodyn,FLBRP-method
#' @examples
#' \dontrun{
#' bd=sim()
#' cpue=stock(bd)
#' setParams(bd) <-cpue
#' setParams(bd)}
#'
setGeneric('setParams<-', function(object,value,...) standardGeneric('setParams<-'))
setMethod('setParams<-', signature(object='biodyn',value='data.frame'), function(object,value) {
nms=c(biodyn:::modelParams(as.character(object@model)),'b0')
object@params=object@params[nms]
object@params =setQ(object,value)
return(object)})
setMethod('setParams<-', signature(object='biodyn',value='FLPar'), function(object,value) {
object@params=value
return(object)})
setMethod('setParams<-', signature(object='biodyn',value='FLQuant'), function(object,value) {
nms=c(biodyn:::modelParams(tolower(as.character(object@model))),'b0')
#value=FLCore::apply(value,2,mean)
object@params =setQ(object,value)
return(object)})
setMethod('setParams<-', signature(object='biodyn',value='FLQuants'), function(object,value,msy=TRUE) {
nms=c(biodyn:::modelParams(as.character(object@model)),'b0')
object@params=object@params[nms]
object@params =setQ(FLCore::iter(object,1),FLQuants(lapply(value, function(x) FLCore::apply(x,2,mean,na.rm=T))))
nms=dimnames(params(object))$param
n = as.numeric(summary(substr(nms,1,1)=='q')['TRUE'])
nms[substr(nms,1,1)=='q'] =paste('q', seq(n),sep='')
nms[substr(nms,1,5)=='sigma']=paste('sigma',seq(n),sep='')
dimnames(params(object))$params=nms
return(object)})
getP<-function(bmsy,k,p=c(0.001,5)){
optimise(function(p,bmsy,k)
(bmsy-k*(1/(1+p))^(1/p))^2, p, bmsy=bmsy,k=k)$minimum}
getRK<-function(msy,bmsy,k,p=c(0.0001,5)){
p=getP(bmsy,k,p)
k=bmsy/((1/(1+p))^(1/p))
r=msy/(k*(1/(1+p))^(1/p+1))
FLPar(c(r=r,k=k,p=p))}
#getRK(100,1000,2000)
setMethod('setParams<-', signature(object='biodyn',value='FLBRP'), function(object,value,msy=TRUE) {
#if (!(is(value)%in%'FLBRP')) return(object)
if (msy){
params(object)[c("r","k","p")]=getRK(value@refpts["msy", "yield"],
value@refpts["msy", "biomass"],
value@refpts["virgin","biomass"])
object@priors[c("r","k","p"),"a"]=params(object)[c("r","k","p")]
}else{
require(popbio)
params(object)["k"]=value@refpts["virgin","biomass"]
params(object)["p"]=getP(c(value@refpts["msy","biomass"]),
c(params(object)["k"]),c(0.001,5))
params(object)["r"]=log(lambda(leslie(eql,fbar=c(value@refpts["crash","harvest"]))))
}
object=fwd(object,catch=catch(object)[,-1])
return(object)})
#' setControl<-
#'
#' @description Sets the control slot in a biodyn object given the parameters in the \code{params}
#' slot. The starting values \code{val} are set to those in \code{params} and the \code{min} and
#' \code{max} bounds to .1 and 10 times these.
#'
#' @param object an object of class \code{biodyn}
#' @param value \code{params} object
#' @param ... any other parameter
#'
#' @seealso \code{\link{controlFn}}
#'
#' @export
#' @rdname setControl
#'
#' @aliases setControl<-,biodyn,FLPar-method setControl<-,biodyn,FLQuant-method setControl<-,biodyn,FLQuants-method
#'
#' @examples
#' \dontrun{
#' bd=sim()
#' setControl(bd)=params(bd)
#' }
#'
setGeneric('setControl<-', function(object,...,value) standardGeneric('setControl<-'))
setMethod('setControl<-', signature(object='biodyn',value='FLPar'), function(object,value,min=0.1,max=10.0) {
if (dims(value)$iter>1 & dims(object@control)$iter==1)
object@control=propagate(control(object),dims(value)$iter)
ctr=object@control
nms=dimnames(object@params)$params
object@control=FLPar(array(rep(c(1,NA,NA,NA),each=length(nms)), dim=c(length(nms),4,dims(value)$iter), dimnames=list(params=nms,option=c('phase','min','val','max'),iter=seq(dims(value)$iter))))
nms.=nms[nms %in% dimnames(ctr)$params]
nits=seq(dims(value)$iter)
object@control[nms.,'phase',nits]=ctr[nms.,"phase",nits]
nms.=nms[nms %in% dimnames(object@params)$params]
object@control[,"val",nits]=object@params[nms.,nits]
object@control[,"min",nits]=object@params[nms.,nits]*min
object@control[,"max",nits]=object@params[nms.,nits]*max
if (!is.na(any(value[nms]<0)) & any(value[nms]<0))
object@control[nms[value[nms]<0],c('min','max')]=object@control[nms[value[nms]<0],c('max','min')]
prr=object@priors
nms=c(nms[!substr(dimnames(object@params)$params,1,1) %in% c('s','q')],
'msy','bmsy','fmsy',
nms[substr(dimnames(object@params)$params,1,1) %in% c('s','q')])
object@priors=array(rep(c(0,0,0.3,1),each=length(nms)), dim=c(length(nms),4), dimnames=list(params=nms,c('weight','a','b','type')))
nms=dimnames(prr)$params[dimnames(prr)$params %in% dimnames(object@priors)$params]
object@priors[nms,]=prr[nms,]
return(object)})
.calcSigma <- function(obs,hat=rep(0,length(obs)),na.rm=T){
n =length(obs[!is.na(obs+hat)])
SS =sum((obs-hat)^2,na.rm=na.rm)
return((SS/n)^.5)}
calcB0<-function(index,q,k,nyrB0=3,error='log'){
if (is.null(nyrB0)) return(params['b0'])
if (error=='log'){
t.<-sweep(log(index[,1:nyrB0,,,,,drop=FALSE]),c(1,6),q,'/')
return(qmax(qmin(exp(apply(t.,c(1,6),mean))/k,1),0))}
if (error=='normal'){
t.<-sweep(index[,1:nyrB0,,,,,drop=FALSE],c(1,6),q,'/')
return(qmax(qmin(apply(t.,c(1,6),mean)/k,1),0))}
}
calcQ<-function(stock,index,error='log',na.rm=T){
stock<-(stock[-length(stock)]+stock[-1])/2
n <-length(stock)
index<-index[seq(n)]
if (na.rm)
n=length(seq(n)[!is.na(index+stock)])
res=switch(error,
normal={q =sum(stock*index, na.rm=T)/sum(stock*stock, na.rm=na.rm)
sigma=.calcSigma(index/(q*stock))
data.frame(q=q,sigma=sigma)
},
log ={q =exp(sum(log(index)-log(stock), na.rm=na.rm)/n)
sigma=.calcSigma(log(index),log(q*stock))
data.frame(q=q,sigma=sigma)},
cv ={res <-sum(index/stock)
sigma2<-.calcSigma(res,na.rm=na.rm)
q <-(-res+(res^2+4*length(index)*sigma2*sum((index/stock)^2)))/(2*length(index)*sigma2)
data.frame(q=q,sigma=sigma)})
return(res)}
setQ=function(object,value,error='log'){
fn=function(value,stock){
if (dims(value)$iter==1)
dimnames(value)$iter=1
if (dims(stock)$iter==1)
dimnames(stock)$iter=1
if (dims(stock)$iter==1 & dims(value)$iter>1)
stock=propagate(stock,dims(value)$iter)
model.frame(mcf(FLQuants(stock=stock,value=value)))}
res=switch(is(value)[1],
FLQuant ={res=fn(value,stock(object));data.frame(name=1,res)},
#FLQuants =ldply(value, fn, model.frame(mcf(FLQuants(stock=stock,value=x))),stock=stock(object)),
FLQuants =ldply(value, function(x,stock) fn(x,stock), stock=stock(object)),
data.frame=merge(model.frame(FLQuants('stock'=stock(object))),value,by='year',all=T))
res=res[!is.na(res$iter),]
if (!('name' %in% names(res)))
names(res)[1]='name'
res=res[!is.na(res$name),]
res=ddply(res, .(name,iter), function(x,log) data.frame(biodyn:::calcQ(x$stock,x$value)),log='log')
its=max(as.numeric(ac(res$iter)))
res.=transform(melt(res,id=c('name','iter')),params=paste(variable,name,sep=''))[,c('params','value','iter')]
names(res.)[2]='data'
#bug
res=as(res.,'FLPar')[,1]
#res=FLCore::iter(res,seq(its))
units(res)='NA'
res=res.[with(res.,order(iter,params)),]
if (dims(object@params)$iter==1)
object@params=propagate(object@params,its)
t.=rbind(object@params,FLPar(as.FLQuant(res)[,1,drop=T]))
dmns=dimnames(t.)
names(dmns)=c('params','iter')
t.=FLPar(array(t.,dim=unlist(lapply(dmns,length)),dimnames=dmns))
units(t.)='NA'
object@params=t.
#object@params=FLPar(rbind(FLPar(object@params),FLPar(res)))
object@params}
setMethod('setControl<-', signature(object='biodyn',value='FLQuant'), function(object,value,
min=0.1,max=10.0) {
setParams(object)<-value
setControl(object,min=min,max=max)<-params(object)
return(object)})
setMethod('setControl<-', signature(object='biodyn',value='FLQuants'),
function(object,value,min=0.1,max=10.0) {
setParams(object)<-value
setControl(object,min=min,max=max)<-params(object)
return(object)})
setMethod('control<-', signature(object='biodyn',value='FLPar'), function(object,value) {
object@control=value
return(object)})
setMethod('control', signature(object='biodyn'), function(object) {
object@control})
#' controlFn
#' @description A utility function to help set up the \code{control} slot in \code{biodyn}
#'
#' @param r a \code{numeric} value with best guess
#' @param k a \code{numeric} value with best guess
#' @param p a \code{numeric} value with best guess, default=1
#' @param b0 a \code{numeric} value with best guess default=1
#' @param phaseR a \code{numeric} value for phase, default=1
#' @param phaseK a \code{numeric} value for phase, default=1,
#' @param phaseP a \code{numeric} value for phase, default=-1,
#' @param phaseB0 a \code{numeric} value for phase, default=-1,
#' @param min a \code{numeric} a multipler for the best guess
#' @param max \code{numeric} a multipler for the best guess
#'
#' @return a \code{control} object
#'
#' @export
#' @rdname controlFn
#'
#' @seealso \code{\link{biodyn}} \code{\link{control}}
#'
#' @examples
#' \dontrun{
#' sim()
#' }
## utility function for setting control object
controlFn=function(r, k, p=1, b0=1,
phaseR=1,phaseK=1,phaseP=-1,phaseB0=-1,
min=.5, max=2){
dmns=list(params=c('r','k','p','b0'),
option=c('phase','min','val','max'),
iter =1)
res=FLPar(array(0,unlist(laply(dmns,length)),dmns))
res[,'val'][]=c(r,k,p,b0)
res[,'min']=res[,'val']*min
res[,'max']=res[,'val']*max
res[,'phase']=c(phaseR,phaseK,phaseP,phaseB0)
res}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.