R/aspic-profile.R

In laurieKell/mpb: Management Procedures

utils::globalVariables('rbind.fill')

#' @title profile
#'
#' @description 
#' Performs a profile using residual sum of squares, fixes some parameters for a range of values 
#' and then estimate the others 
#'
#' @param fitted: an \code{aspic} object
#' @param which: \code{character} giving the parameters to do the profile for, i.e. to fix.
#' @param range; \code{numeric} relative values by which to vary parameter values, default seq(0.5,1.5,length.out=21). 
#' @param fn: \code{function} that gives values to be profiled.
#' @param run: \code{logical} if \code{TRUE} then returns profile, otherwise it just sets the control object-
#' 
#' @return a \code{data frame} with results turned by \code{fn} by values in \code{which}. 
#' @seealso \code{\link{biodyn},\link{fit}}
#'
#' @export
#' @docType methods
#' @rdname profile
#'
profileFn=function(fitted,which,
                   range=seq(0.5,1.5,length.out=21),
                   fn   =function(x) cbind(model.frame(params(x)), 
                                           model.frame(refpts(x)),
                                           model.frame(x@objFn)[,-3]),
                   #ll=x@ll[,2,drop=T]/(x@ll[,3,drop=T]-1)),
                   run=TRUE,
                   min=0.8,max=1.2,...){
  
  if (dims(fitted)$iter>1) stop("can only be done for a single iter")
  
  if (length(range)==1) range=c(range,2-range)
  
  if (dim(fitted@control)[3]==1){
    sq=list(range)
    sq=do.call("expand.grid",sq[rep(1,length(which))])
    names(sq)=which
    
    setControl(fitted,min=min,max=max)=params(fitted)
    fitted@control=propagate(fitted@control,dim(sq)[1])
    
    #fix bounds for msy & k
    for (i in which){
      fitted@control[i,"val"]=params(fitted)[i,]*sq[,i]
      fitted@control[i,"min"]=min(fitted@control[i,"val"])*0.9
      fitted@control[i,"max"]=max(fitted@control[i,"val"])*1.1
      fitted@control[i,"fit"]=0
      
      if (("k"  == i) & all(c("msy","k")%in%dimnames(fitted@control)$params))
        fitted@control["msy","min"]=min(fitted@control["k","min"],fitted@control["msy","min"])*0.9
      if (("msy"== i) & all(c("msy","k")%in%dimnames(fitted@control)$params))
        fitted@control["msy","min"]=min(fitted@control["k","min"],fitted@control["msy","min"])*0.9
    }
    
    if (!run) return(fitted)
    
    fitted=fit(fitted)
    fitted=fn(fitted)}
  else
    fitted@control=profileGrid(fitted@control,which,range)
  
  return(fitted)}

setMethod("profile", signature(fitted="aspic"), 
          function(fitted,which,
                   range=seq(0.5,1.5,length.out=21),
                   fn   =function(x) cbind(model.frame(params(x)), 
                                           model.frame(refpts(x)),
                                           model.frame(x@objFn)[,-3]),
                                           #ll=x@ll[,2,drop=T]/(x@ll[,3,drop=T]-1)),
                   run  =TRUE,min=0.8,max=1.2,...)
     profileFn(fitted=fitted,which=which,range=range,fn=fn,run=run,min=min,max=max))

setMethod('profile',  signature(fitted='aspics'),
          function(fitted,which,
                   range=seq(0.5,1.5,length.out=21),
                   fn   =function(x) cbind(model.frame(params(x)), 
                                           model.frame(refpts(x)),
                                           model.frame(x@objFn)[,-3],
                                           t(unlist(x@ll)[1,drop=T])),
                   run=TRUE,
                   .multicombine=T,.maxcombine=10,.packages=c("aspic","plyr")){
            
            if (run) .combine=rbind.fill else .combine=list

            res=foreach(i=names(fitted), .combine=.combine,
                        .multicombine=.multicombine,
                        .maxcombine  =.maxcombine,
                        .packages    =.packages) %dopar% {  
                        
                         profile(fitted[[i]],which=which,range=range,fn=fn,run=run)}
            
            if (!run) {
              res       =aspics(res)
              names(res)=names(fitted)}
            else {
              ns =ldply(fitted, function(x) length(unique(index(x)[,"name"])))[,2]
              ns =ns*length(range)*length(which)
              ns =rep(seq(length(ns)),ns)
              res=cbind(.id=names(fitted)[ns],res)
              }
            
            res})

#' @title fnProfile
#' 
#' @description 
#' Outputs as a data.frame a summary of parameters and RSS etc by data component 
#' generated when doing a profile
#' 
#' @param x: an \code{aspic} object
#' 
#' @return a \code{data frame} with results by data component. 
#' @seealso \code{\link{biodyn},\link{profile} \link{fit}}
#'
#' @export
#' @docType methods
#' @rdname fnProfile 
#'
#' @examples
#' \dontrun{
#' data(asp)
#' dcK=profile(asp,which=c("k"),range=seq(0.2,2.0,length.out=21),fn=fnProfile)
#' ggplot(dcK)+geom_line(aes(k,value,group=variable,col=variable))+
#'  theme_ms(12,legend.position="bottom")+
#'  ylab("Residual Sum of Squares")+xlab("K")
#' 
#' }       
fnProfile=function(x) {
  
  res=cbind(model.frame(params(x),drop=T)[,-(dim(params(x))[1]+1)],
            model.frame(refpts(x),drop=T)[,-c(1,4)],
            model.frame(FLQuants(stock  =stock(  x)[,ac(range(x)["maxyear"])]%/%bmsy(x),
                                 harvest=harvest(x)[,ac(range(x)["maxyear"])]%/%fmsy(x)),drop=T)[,-1],
            model.frame(x@ll[,"ss"])[,seq(dim(x@ll)[1])]
  )
  
  nms=dimnames(params(x))$params
  res=melt(res,id=c(dimnames(refpts(x))$refpts[-1],nms,c("stock","harvest")))
  
  res=transform(res,name=unique(index(x)$name)[res$variable])
  res}


#levels=(logLik(nsher)-qchisq(.95,2)/2), col="navy")

# ### debugging stuff
# data(bd)
# fitted=biodyn(factor("pellat"),params(bd),catch=catch(bd))
# cpue=rlnorm(1,log(stock(bd)),.2)[,-60]
# setParams(fitted)     =cpue
# 
# 
# attach(list(length.out=11, range=0.5, ci=c(0.25, 0.5, 0.75, 0.95),
#             plot=TRUE,fixed=c()))
# which="r"
# fixed=c("p","b0")
# ###
# rtn=profile(swon[[1]],which="k",fixed="b0",length.out=31,range=c(.75,1.5))
# ggplot(rtn)+geom_line(aes(msy,rss))

# \dontrun{
#  data(asp)
#  res=profile(asp,which="msy",range=seq(0.75,1.25,length.out=21))
#  ggplot(res)+geom_line(aes(k,rss))
# 
#  fn =function(x) cbind( data.frame(model.frame(params(x)),rss=sum(diags(x)$residual^2,na.rm=T)))
#  res=profile(asp,which="msy",length.out=3,range=c(0.5,1.1),fn=fn)
# 
#   fn=function(x) ddply(x@diags, .(name), with, sum(residual^2,na.rm=T)/sum(count(!is.na(residual))))
# }       

# refpts<-function(object){
#   msy=fmsy(object)*bmsy(object)
#   dimnames(msy)$params="msy"
#   rbind(msy,
#         bmsy(object),
#         fmsy(object))}