biodyn: Biomass Dynamic Model

utils::globalVariables(c('admbCor','admbProfile',
                         'daply','residual', 'count'))

utils::globalVariables(c('laply','llply','maply','mlply'))
utils::globalVariables('alply')

diagsFn=function(res){
  res$residualLag <- c(res$residual[-1],NA)
  
  qqLine <- function(x,y){ 
    qtlx <- quantile(x, prob=c(0.25,0.75), na.rm=T)
    qtly <- quantile(y, prob=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)}
  
  try({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}

#### ADMB ###################################################################################
## 1) setExe copies exe from bin and creates a temo dir
## 2) Write data files
## 3) Run exe
## 4) Read output files
#############################################################################################

#tmp=model.frame(FLPar(array(NA,dim=c(3,2,1),dimnames=list(param=c('a','b','c'),var=c('x','y'),iter=1))))
#tmp=as.data.frame(FLPar(array(NA,dim=c(3,2,1),dimnames=list(param=c('a','b','c'),var=c('x','y'),iter=1))))

## copies exe into temp dir ready to run
setExe=function(exeNm,package,dir=tempdir()){
  ##### set up temp dir with exe for data files
  
  # Linux
  if (R.version$os=='linux-gnu') {
    exe = paste(system.file('bin', 'linux', package=package, mustWork=TRUE),exeNm, sep='/')
    if (length(grep("-rwxrwxr-x",system(paste("ls -l",exe),inter=TRUE)))==0)
      warning("Executable privilege not set for \n",exe,call.=FALSE)
    
    file.copy(exe, dir)
    dir = paste(dir, '/', sep='')
     
    # Windows
  } else if (.Platform$OS.type == 'windows') {
    exe = paste(system.file('bin', 'windows', package=package, mustWork=TRUE), paste(exeNm, '.exe', sep=''), sep='/')
    file.copy(exe, dir)
    print(exe)
    print(dir)
    dir = paste(dir, '\\', sep='')
    
    # Mac OSX
  }else 
    stop()
  
  oldwd = getwd()
  
  # change wd to avoid exe case bug
  setwd(dir)
  
  oldwd}

setPella=function(obj, exeNm='pella', dir=tempdir()) {
  # create input files ################################
  dgts=options()$digits
  options(digits=22)

  # cpue
  if (is.FLQuant(obj[[2]])){  
     idx=model.frame(FLQuants(index=obj[[2]]), drop=TRUE)
     idx=data.frame(idx,name=1)
  }else if ('FLQuants' %in% class(obj[[2]])){  
     idx=as.data.frame(obj[[2]], drop=TRUE)
     names(idx)[2:3]=c('index','name')
     idx=transform(idx,name=as.numeric(name))
  }
  idx=idx[!is.na(idx$index),]

  bd.=obj[[1]]
 
  nms=c(biodyn:::modelParams('pellat'),'b0')
  if (length(unique(idx$name))>0)
    nmIdx=paste(c('q','sigma'), rep(unique(idx$name),each=2),sep='')
  else  
    nmIdx=c('q','sigma')

  ctl = bd.@control[nms,]
  ctl[,2:4] = ctl[,c(2,4,3)]
  ctl=alply(array(c(ctl),dim=dim(ctl),dimnames=dimnames(ctl))[,,1,drop=T],1)
#  ctl = alply(ctl,1)
  names(ctl) = nms
#   
#   # ctl file
#   ctl        =bd.@control[nms,,1]
#   ctl[,2:4,1]=ctl[,c(2,4,3),1]
# print("set 1")
#   ctl        =alply(ctl,1)
# print("set 2")
#   names(ctl) =nms

  biodyn:::writeADMB(ctl, paste(dir, '/', exeNm, '.ctl', sep=''),FALSE)

  cat('# q ####################\n', file=paste(dir, '/', exeNm, '.ctl', sep=''),append=TRUE)

  ctl           = bd.@control[nmIdx[grep('q',nmIdx)],,1]
  ctl[,2:4,1]   = ctl[,c(2,4,3),1]
  ctl           = alply(t(matrix(ctl,dim(ctl))),1)
  names(ctl)    = c('phase','lower','upper','guess')

  biodyn:::writeADMB(ctl, paste(dir, '/', exeNm, '.ctl', sep=''),TRUE)
  
  cat('# sigma ################\n', file=paste(dir, '/', exeNm, '.ctl', sep=''),append=TRUE)
  ctl           = bd.@control[nmIdx[grep('s',nmIdx)],,1]
  ctl[,2:4,1]     = ctl[,c(2,4,3),1]
  ctl           = alply(t(matrix(ctl,dim(ctl))),1)
  names(ctl)    = c('phase','lower','upper','guess')

  biodyn:::writeADMB(ctl, paste(dir, '/', exeNm, '.ctl', sep=''),TRUE)
  
  # prr file
  prr = bd.@priors[c(nms,c('msy','bmsy','fmsy'),nmIdx),] 
  prr = alply(prr,1)

  names(prr) = dimnames(bd.@priors)$params
  biodyn:::writeADMB(prr, paste(dir, '/', exeNm, '.prr', sep=''))

  # ref file
  if (is.na(bd.@ref["refyr"])) 
       bd.@ref["refyr"]=as.integer(range(bd.)["maxyear"]-(range(bd.)["minyear"]+range(bd.)["maxyear"])/2)

  biodyn:::writeADMB(bd.@ref, paste(dir, '/', exeNm, '.ref', sep=''))

  # write data
  ctc = as.list(model.frame(FLQuants("catch"=bd.@catch), drop=TRUE))
  ctc = c(nYrs=length(ctc[[1]]), ctc)
  res = c(ctc, c(nIdxYrs=dim(idx)[1], nIdx=length(unique(idx$name)), idx))

  biodyn:::writeADMB(res, paste(dir, '/', exeNm, '.dat', sep=''))

#   # propagate as required
#   its = dims(bd)$iter
#   
#   # params
#   bd@params = propagate(params(bd)[,1], its)
#   # stock
#   stock(bd)  = FLQuant(dimnames=dimnames(stock(bd))[1:5], iter=its)
#   

  # vcov
  vcov(bd.)=FLPar(array(NA, dim     =c(dim(params(bd.))[1],dim(params(bd.))[1],1), 
                           dimnames=list(params=dimnames(params(bd.))[[1]],
                                         params=dimnames(params(bd.))[[1]],iter=1)))
  
  options(digits=dgts)
  return(bd.)}

getPella=function(obj, exeNm='pella') {
  t1 = read.table(paste(exeNm,'.rep',sep=''),skip =18,header=T) 
  # params
  t2 = unlist(c(read.table(paste(exeNm,'.rep',sep=''),nrows=4)))
  q. = unlist(c(read.table(paste(exeNm,'.rep',sep=''),nrows=1,skip=8)))

  s. = unlist(c(read.table(paste(exeNm,'.rep',sep=''),nrows=1,skip=10)))

  nms=c('r','k','b0','p')
  obj@params[nms,] = t2

  obj@params[grep('q',dimnames(obj@params)$params),]=q. 
  obj@params[grep('s',dimnames(obj@params)$params),]=s. 
  t3 = unlist(c(read.table(paste(exeNm,'.rep',sep=''),skip=dim(params(obj))[1]*2,nrows=2,header=F)))

  obj@objFn['ll'] =t3[length(t3)]
  obj@objFn['rss']=t3[length(t3)-1]

  us=paste('u',seq(length(dimnames(params(obj))$params[grep('q',dimnames(params(obj))$params)])),sep='')
  vals=FLPar(unlist(biodyn:::readADMB('lls.txt')),dimnames=list(params=us,iter=1))
  obj@ll=vals

  # stock biomass
  obj@stock[,1:dim(t1)[1]] = unlist(c(t1['stock'])) 

  return(obj)} 

#FLParBug sim@control['r','val',1]=c(.5,.6)
#exe(object, exeNm='biodyn', dir=tempdir(), set=biodyn:::set, get=biodyn::set, cmdOps=paste('-maxfn 500'))
  
activeParams=function(obj) dimnames(obj@control)$params[c(obj@control[,'phase']>-1)]

## runs exe

#' fit
#'
#' Estimates parameters \code{biodyn} class by fitting catch and CPUE indices
#' 
#' @param   object an object of class \code{biodyn}
#' @param   index an \code{FLQuant}, \code{FLQuants} or  \code{data.frame} object with CPUE indices
#' @param   ... other arguments
#'
#' @export
#' @rdname fit
#'
#' @aliases fit fit-method fit,biodyn,FLQuant-method   fit,biodyn,FLQuantJK-method  fit,biodyn,FLQuants-method
#' 
#' @examples
#' \dontrun{
#'    #simulate an object with known properties
#'    bd=sim()
#'    bd=window(bd,end=49)
#'    
#'    #simulate a proxy for stock abundance
#'    cpue=(stock(bd)[,-dims(bd)$year]+stock(bd)[,-1])/2
#'    cpue=rlnorm(1,log(cpue),.2)
#'    
#'    #set parameters
#'    setParams(bd) =cpue
#'    setControl(bd)=params(bd)
#'    control(bd)[3:4,"phase"]=-1
#'    
#'    #fit
#'    bd=fit(bd,cpue)
#' }
setGeneric('fit',   function(object,index,...)     standardGeneric('fit'))
setMethod('fit',signature(object='biodyn',index='FLQuant'),
          function(object,index=index,exeNm='pella',package='biodyn', 
                   dir=tempdir(),
                   set=setPella,
                   get=getPella,cmdOps=paste('-maxfn 500 -iprint 0'))
   fitPella(object,index=index,exeNm=exeNm,package=package, 
            dir=dir,
            set=set,
            get=get,cmdOps=cmdOps))

setMethod('fit',signature(object='biodyn',index='FLQuants'),
          function(object,index=index,exeNm='pella',package='biodyn', 
                   dir=tempdir(),
                   set=setPella,
                   get=getPella,cmdOps=paste('-maxfn 500 -iprint 0'))
            fitPella(object,index,exeNm,package, 
                     dir=dir,
                     set=set,
                     get=get,cmdOps=cmdOps))

setMethod('fit',signature(object='biodyn',index='FLQuantJK'),
          function(object,index=index,exeNm='pella',package='biodyn', 
                   dir=tempdir(),
                   set=setPella,
                   get=getPella,cmdOps=paste('-maxfn 500 -iprint 0')){

            object=propagate(object,dims(index)$iter)
            
            index =as.FLQuant(index)

            object=fitPella(object,index=index,exeNm=exeNm,package=package, 
                     dir=dir,
                     set=set,
                     get=get,cmdOps=cmdOps)
            
            attributes(object)['jk']=TRUE
            
            object})


fitPella=function(object,index=index,exeNm='pella',package='biodyn', 
                  dir=tempdir(),
                  set=setPella,
                  get=getPella,cmdOps=paste('-maxfn 500 -iprint 0'))          
  {
  first=TRUE   
  
  catch=NULL
  if (dims(object)$iter==1 & 1<ifelse('FLQuant'%in%is(index),dims(index)$iter>1,max(laply(index,function(x) dims(x)$iter))))
     {
     catch=catch(object)

    if ("FLQuants"%in% is(index))
      its=max(laply(u,function(x) dims(x)$iter))
    else
      its=dims(index)$iter
     catch(object)=propagate(catch(object),its)   
    }

  max=min(dims(catch(object))$maxyear,ifelse('FLQuant'%in%is(index),dims(index)$maxyear,max(laply(index,function(x) dims(x)$maxyear))))
  if (!is.na(range(object)['maxyear'])) max=min(max,range(object)['maxyear']) 
  min=min(dims(catch(object))$minyear,ifelse('FLQuant'%in%is(index),dims(index)$minyear,max(laply(index,function(x) dims(x)$minyear))))
  if (!is.na(range(object)['minyear'])) min=max(min,range(object)['minyear'])
  object=window(object,start=min,end=max)

  if ('FLQuant'%in%is(index)){ 
    index =window(index,start=min,end=max)
  }else if ('FLQuants'%in%is(index)){
    index=FLQuants(llply(index, window,start=min,end=max))}
  
  slts=getSlots('biodyn')
  slts=slts[slts %in% c('FLPar','FLQuant')]

  #oldwd =setExe(exeNm,package,dir)
  oldwd=getwd()
  setwd(dir)
  exe()

  object=list(object,index)
  bd =object[[1]]
  its=max(maply(names(slts), function(x) { 
          #print(x)
          #print(dims(slot(bd,x))$iter)
          dims(slot(bd,x))$iter 
          })) 
  its=max(its,dims(bd@control)$iter)

  nms=dimnames(params(bd))$params
  bd@vcov   =FLPar(array(as.numeric(NA), dim=c(length(nms),length(nms),its), dimnames=list(params=nms,params=nms,iter=seq(its))))
  bd@hessian=bd@vcov
   
  us=paste('u',seq(length(dimnames(params(bd))$params[grep('q',dimnames(params(bd))$params)])),sep='')
  bd@ll=FLPar(NA,dimnames=list(params=us,iter=seq(1)))

  if (its>1){
   
      ## these are all results, so doesnt loose anything
      bd@stock  =FLCore::iter(bd@stock,  1)
      bd@params =FLCore::iter(bd@params, 1)
      bd@objFn  =FLCore::iter(bd@objFn,  1)
      bd@vcov   =FLCore::iter(bd@vcov,   1)
      bd@ll     =FLCore::iter(bd@ll,     1)
      bd@hessian=FLCore::iter(bd@hessian,1)
      bd@mng    =FLPar(a=1)
      bd@mngVcov=FLPar(a=1,a=1)
      
      if (dim(bd@stock)[6]==1) bd@stock=propagate(bd@stock, iter=its, fill.iter=TRUE)      
      if (dim(bd@catch)[6]==1) bd@stock=propagate(bd@catch, iter=its, fill.iter=TRUE)     
      
      pnms   <- getSlots(class(bd))
      pnames <- names(pnms)[pnms == 'FLPar']
      for(i in pnames){
        slot(bd, i)=FLCore::iter(slot(bd, i),1)
        slot(bd, i) <- propagate(slot(bd, i), its)}
      
      #bd=propagate(bd,its)      
      }
   
  #print(slot(object[[1]],'control'))

  cpue=object[[2]]
  bd2 =object[[1]]

  for (i in seq(its)){       
     object[[2]] = FLCore::iter(cpue,i) 
     
     for (s in names(slts)[-(7:8)]){      
        slot(object[[1]],s) = FLCore::iter(slot(bd2,s),i) 
        }  
     
     object[[1]]=set(object,exeNm,dir)

     exe = paste(system.file('bin', 'linux', package="biodyn", mustWork=TRUE),exeNm, sep='/')      

    #bug in windows
    try(
      if (length(grep("-rwxrwxr-x",system(paste("ls -l",exe),inter=TRUE)))==0)
        warning("Executable privilege not set for \n",exe,call.=FALSE) )

     # run
     #system(paste('./', exeNm, ' ', cmdOps, sep=''))
     system(paste(exeNm, ' ', cmdOps, sep=''))
    
     # gets results
     object[[1]]=getPella(object[[1]], exeNm)     
     s=names(slts)[slts=='FLQuant']
     for (s in s[!(s=="catch")])
       try(FLCore::iter(slot(bd,s),i) <- slot(object[[1]],s))
          
     if (its<=1 & file.exists(paste(dir,'admodel.hes',sep='/'))){
       ##hessian
       x<-file(paste(dir,'admodel.hes',sep='/'),'rb')
       nopar<-readBin(x,'integer',1)
       H<-matrix(readBin(x,'numeric',nopar*nopar),nopar)
       try(bd@hessian@.Data[activeParams(object[[1]]),activeParams(object[[1]]),i] <- H, silent=TRUE)
       close(x)
     
       ## vcov
       #print(file.exists(paste(dir,'admodel.cov',sep='/')))
       if (file.exists(paste(dir,'admodel.cov',sep='/')))
         try(bd@vcov@.Data[activeParams(object[[1]]),activeParams(object[[1]]),i] <- cv(paste(dir,'admodel.hes',sep='/')), silent=TRUE) 
        
       #if (file.exists(paste(dir,'admodel.cov',sep='/'))){
       #   x<-file(paste(dir,'admodel.cov',sep='/'),'rb')
       #   nopar<-readBin(x,'integer',1)
       #   H<-matrix(readBin(x,'numeric',nopar*nopar),nopar)
       #   try(bd@vcov@.Data[activeParams(object[[1]]),activeParams(object[[1]]),i] <- H, silent=TRUE)
       #close(x)}
       
       if (file.exists(paste(dir,'pella.hst',sep='/')))
          bd@hst=biodyn:::admbProfile(paste(dir,'pella.hst',sep='/'))$profile
       if (file.exists(paste(dir,'lpr.plt',sep='/')))
          bd@profile=mdply(data.frame(var=c("r", "k","bnow","fnow",
                                            "bnow","fnow","bnowthen","fnowthen",
                                            "msy","bmsy","fmsy","cmsy",
                                            "bmsy","ffmsy","bk","fr",
                                            "bratio","fratio","slopeb","slopef")),
                                function(var){
                                    #print(var)
                                    fl=paste("lp",var,".plt",sep="")
                                    if (file.exists(fl))
                                      biodyn:::admbPlt(fl)})
       }

     bd@params@.Data[  ,i] = object[[1]]@params
     bd@control@.Data[,,i] = object[[1]]@control
     #bd@objFn@.Data[   ,i] = object[[1]]@objFn

     #FLParBug
     bd@ll@.Data[,i][] = unlist(c(object[[1]]@ll))

     if (file.exists('pella.std')){
       err1=try(mng.<-read.table('pella.std',header=T)[,-1])
    
       err2=try(mngVcov.<-fitFn(paste(dir,'pella',sep='/'))$vcov)

     ## FLPar hack
     if (first) {

       if (any(is(err1)!='try-error')) 
         bd@mng=FLPar(array(unlist(c(mng.[   ,-1])), dim     =c(dim(mng.)[1],2,its),
                                                     dimnames=list(param=mng.[,1],var=c('hat','sd'),iter=seq(its))))
        
       if (any(is(err2)!='try-error')) 
         bd@mngVcov<-FLPar(array(unlist(c(mngVcov.)),dim     =c(dim(mng.)[1],dim(mng.)[1],its),
                                                     dimnames=list(param=dimnames(mngVcov.)[[1]],
                                                                   param=dimnames(mngVcov.)[[1]],iter=seq(its))))
 
       first=!first  
    }else{     
       try(if (is(err1)!='try-error') bd@mng@.Data[,,i][]=unlist(c(mng.[,-1])))
       try(if (is(err2)!='try-error') bd@mngVcov@.Data[,,i][]=unlist(c(mngVcov.)))
       }}
  }
      
  units(bd@mng)='NA'

  bd=fwd(bd,catch=catch(bd)[,rev(dimnames(catch(bd))$year)[1]])

  if (length(grep('-mcmc',cmdOps))>0 & length(grep('-mcsave',cmdOps))>0){
    #'-mcmc 100000 -mcsave 100'
     setMCMC=function(obj,dir){
       ps=read.psv(paste(dir,'pella.psv',sep='/'))

       dmns=list(params=activeParams(obj),iter=seq(dim(ps)[1]))

       ps=array(t(ps),dim=unlist(llply(dmns,length)),dimnames=dmns)
       ps=FLPar(ps)
       
       units(ps)='NA'
       ps}
     
    par=setMCMC(bd,dir)
    
    cmd=strsplit(cmdOps,',')
    grp=unlist(gregexpr('-mcmc',cmd[[1]])) 
    mcmc =sub(' +', '', cmd[[1]][grp>0]) 
    mcmc =as.numeric(substr(mcmc,6,nchar(mcmc)))
    grp=unlist(gregexpr('-mcsave',cmd[[1]])) 
    mcsave=sub(' +', '', cmd[[1]][grp>0])
    mcsave=sub(' +', '', mcsave)
    mcsave=as.numeric(substr(mcsave,8,nchar(mcsave)))
     
    bd@params=propagate(bd@params[,1],dims(par)$iter)
    bd@params[dims(par)$params,]=par
    bd@stock=propagate(bd@stock,dim(params(bd))[2])
    bd=fwd(bd,catch=catch(bd))  

    attributes(bd@params)[['mcmc']]  =mcmc
    attributes(bd@params)[['mcsave']]=mcsave 
    }
 
  if ("FLQuant"%in%class(index)) index=FLQuants("1"=index)

  if (its<=1){
      i=0
      bd@diags=dgs=ldply(index,function(index){
        i<<-i+1
        stockHat=(stock(bd)[,-dims(stock(bd))$year]+stock(bd)[,-1])/2
        hat     =stockHat*params(bd)[paste("q",i,sep="")]
        res=model.frame(mcf(FLQuants(
          index   =index,
          stock   =stock(bd),
          stockHat=stockHat,
          hat     =hat,
          residual=log(index/hat))),drop=T)
        
        diagsFn(res)})}

  setwd(oldwd)

#  if (!is.null(catch)) catch(object)=catch


  return(bd)}
            
#library(matrixcalc)
#is.positive.definite(H)

#vc=vcov(bd)[activeParams(bd),activeParams(bd),drop=T]
#vc[lower.triangle(vc)==0]=t(vc)[lower.triangle(vc)==0]
#is.positive.definite(vc,tol=1e-8)
#(vc-t(vc))/vc


admbCor=function(fl='pella.cor'){
  if (!file.exists(fl)) return
  
  tmp=options()
  options(warn=-1)
  
  res=scan(fl,sep='\n',what=as.character(),quiet=TRUE)[-1]
  res=maply(res, str_trim)
  res=strsplit(res, ' +')
     
  nms=laply(res,function(x) x[2])[-1]
  hat=laply(res,function(x) as.numeric(x[3]))[-1]
  sd =laply(res,function(x) as.numeric(x[4]))[-1]
  
  npar=length(res)-1
  
  val =as.numeric(unlist(llply(res[-1], function(x) strsplit(x,' +')[-(1:4)])))
  cor=upper.triangle(array(1,dim=c(npar,npar),dimnames=list(nms,nms)))
  cor[upper.triangle(cor)==1]=val
  cor=cor+t(cor)
  diag(cor)=1
  
  vcov=cor*sd%o%sd
  
  options(tmp)
  
  hat =FLPar(hat, units='NA')
  dimnames(hat)$params=dimnames(vcov)[[1]]
  vcov=FLPar(vcov,units='NA')
  names(vcov)[1:2]=c('params','params')
  
  
  return(FLPars(hat =hat,
                vcov=vcov))}

getDiags=function(fl='pella.rep'){
  skip=grep('Model',scan(fl,sep='\n',what=as.character(),quiet=TRUE))
      
  res=read.table(fl,skip=skip,header=TRUE)
  res=transform(res,residual=log(index/hat))
  res=diagsFn(res)
  
  res$harvest=res$catch/res$stock
  res$stock  =res$stock.
   
  res[,-7]}

calcElasticity=function(bd,mn=3,rg=5){
  
  elasFn=function(x,dmns,bd,mn,rg) {
    
    params(bd)[dmns]=exp(x)
    bd=fwd(bd,catch=catch(bd))
    
    maxYr =ac(range(bd)['maxyear'])
    max1  =ac(range(bd)['maxyear']-(seq(mn)-1))
    max2  =ac(range(bd)['maxyear']-(seq(mn)-1+mn))
    maxR  =ac(range(bd)['maxyear']-(seq(rg)-1))
     
    smy=c(  stock(bd)[,maxYr],
          harvest(bd)[,maxYr],
              msy(bd),
             fmsy(bd),
             bmsy(bd),
            catch(bd)[,maxYr]/msy( bd),
            stock(bd)[,maxYr]/bmsy(bd),
            harvest(bd)[,maxYr]/fmsy(bd),
            mean(stock(bd)[,max1])/mean(stock(bd)[,max2]),
            mean(harvest(bd)[,max1])/mean(harvest(bd)[,max2]),
            abs(coefficients(lm(data~year,as.data.frame(stock(  bd)[,maxR],drop=T)))['year']),
            abs(coefficients(lm(data~year,as.data.frame(harvest(bd)[,maxR],drop=T)))['year']))     
      
    return(log(smy))}
  
  parNms=c(biodyn:::modelParams(model(bd)),'b0')

  jbn=jacobian(elasFn,log(c(bd@params[parNms])),dmns=parNms,bd=bd,mn=mn,rg=rg)
  
  
  dimnms=list(params=parNms,
              stat  =c('stock',    'harvest',
                       'msy',      'bmsy',     'fmsy',
                       'catchMsy', 'stockBmsy','harvestFmsy',
                       'stockMn',  'harvestMn',
                       'stockRg',  'harvestRg'))
  
  jbn=FLPar(array(t(jbn),dim=dim(t(jbn)),dimnames=dimnms))
  units(jbn)='NA'
  
  return(jbn)}

exe=function(package='biodyn'){

  sep =  function() if (R.version$os=='linux-gnu') ':' else if (.Platform$OS=='windows') ';' else ','
  
  # Linux
  if (R.version$os=='linux-gnu') {
    path = paste(system.file('bin', 'linux',  package=package, mustWork=TRUE), sep='/')
    # Windows
  } else if (.Platform$OS.type == 'windows') {
    
    path = paste(system.file('bin', 'windows', package=package, mustWork=TRUE), sep='/')
    # Mac OSX
  }else 
    stop()
  
  #path='C:\\R\\R-2.15.1\\library\\aspic\\bin\\windows'
  
  path <- paste(path, sep(), Sys.getenv('PATH'),sep='')
   
  Sys.setenv(PATH=path)

return(path)}
#asp=aspic('http://gbyp-sam.googlecode.com//svn//trunk//tests//aspic//swon//2009//high//aspic.inp')
#asp=fit(asp)

calcSS=function(x) daply(x@diags, .(name), 
                         with, sum(residual^2,na.rm=T)/sum(count(!is.na(residual))))

#FLParBug
#tst=FLPar(as.numeric(NA),dimnames=list(params=c('a','b'),col=c('x','y'),iter=1:2))
#tst['a','x',2]=3

fitFn=function(file){

  res=biodyn:::admbFit(file)

  #est        
  hat =FLPar(array(c(res$est),dim     =c(length(res$names),1),
                              dimnames=list(params=res$names,iter=1)))
  units(hat)='NA'  
    
  #std        
  std =FLPar(array(c(res$std),dim     =c(length(res$names),1),
                              dimnames=list(params=res$names,iter=1)))
  units(std)='NA'  
    
  #cor 
  cor =FLPar(array(c(res$cor),dim     =c(length(res$names),length(res$names),1),
                              dimnames=list(params=res$names,params=res$names,iter=1)))
  units(cor)='NA'  
    
  #cov
  vcov=FLPar(array(c(res$cov),dim     =c(length(res$names),length(res$names),1),
                              dimnames=list(params=res$names,params=res$names,iter=1)))
    
  units(vcov)='NA'  
    
  return(list(std       =std,hat=hat,cor=cor,vcov=vcov,  
              nlogl     =res$nlogl,      
              maxgrad   =res$maxgrad,    
              npar      =res$npar,       
              logDetHess=res$logDetHess))}

#file='/tmp/Rtmp6dBbkc/pella'
#fitFn(file)