R/biodyn-F.R

In laurieKell/biodyn: Biomass Dynamic Model

globalVariables(c("slv","ucminf"))

#PELLA , J. J., 1967 A study of methods to estimate the Schaefer model parameters with special reference to the yellowfin tuna fishery in the eastern tropical Pacific ocean. University of Washington, Seattle.
#PELLA , J. J., and P. K. T OMLINSON , 1969 A generalized stock production model. Bulletin of the Inter-American Tropical Tuna Commission 13: 419-496.
#PRAGER , M. H., 1994 A suite of extensions to a nonequilibrium surplus-production model. U. S. Fishery Bulletin 92: 374-389.

alpha=function(r,F) r-F
beta =function(r,K) r/K

yield<-function(F,B,r,K)
    (F/(r/K))*log(1-(r/K)*B*(1-exp((r-F)))/(r-F))

fnY   =function(F,C,B,r,K)
  C-yield(F,B,r,K)

minY   =function(F,C,B,r,K)
  fnY(F,C,B,r,K)^2

gradY =function(F,C,B,r,K){
  .expr1  <- r/K
  .expr2  <- F/.expr1
  .expr3  <- .expr1*B
  .expr4  <- r-F
  .expr5  <- exp(.expr4)
  .expr7  <- .expr3 * (1 - .expr5)
  .expr9  <- 1 - .expr7/.expr4
  .expr10 <- log(.expr9)
  .value  <- C - .expr2 * .expr10
  -(1/.expr1 * .expr10 - .expr2 * ((.expr3 *.expr5/.expr4 + .expr7/(.expr4^2))/.expr9))
  }

gradMinY=function(F,C,B,r,K)
  2*minY(F,C,B,r,K)*gradY(F,C,B,r,K)

fn<-function(F,C,B,r,K)
  C*(r/K)/(log(r/K*B*(exp(r-F)-1)/(r-F)+1))

fnF   =function(F,C,B,r,K)
  F-fn(F,C,B,r,K)

minF   =function(F,C,B,r,K)
  fnF(F,C,B,r,K)^2

gradF=function(F,C,B,r,K){
  .expr1 <- r/K
  .expr2 <- C*.expr1
  .expr3 <- .expr1*B
  .expr4 <- r-F
  .expr5 <- exp(.expr4)
  .expr7 <- .expr3*(.expr5 - 1)
  .expr9 <- .expr7/.expr4 + 1
  .expr10 <- log(.expr9)
  1 - .expr2 * ((.expr3*.expr5/.expr4 - .expr7/(.expr4^2))/.expr9)/.expr10^2
  }

gradMinF=function(F,C,B,r,K)
  2*minF(F,C,B,r,K)*gradF(F,C,B,r,K)

NewRhap=function(x, func, grad)
  x - func / grad

prj<-function(F,C,B,r,K)
  alpha(r,F)*B*exp(alpha(r,F))/(alpha(r,F)+beta(r,K)*B*(exp(alpha(r,F))-1))

nr<-function(C,B,r,K,B0,tolVal=1e-10,niter=200,yield=TRUE){
  
  F=-log(1-C/B[,-dim(B)[2]])*.5
  for (t in rev(rev(dimnames(C)$year)[-1])){
       iters=0;tol=1;val=1
       
       while (abs(val)>tolVal&iters<niter){
         iters  =iters+1
         if (yield){
           func   = fnY(F[,t],C[,t],B[,t],r,K) 
           grad   = gradY(F[,t],C[,t],B[,t],r,K)
         }else{
           func   = fnF(F[,t],C[,t],B[,t],r,K) 
           grad   = gradF(F[,t],C[,t],B[,t],r,K)
           }

         val=func/grad
        
         F[,t]=F[,t]-val
         
       #print(c(func,val,c(F[,t])))
       }
       
   #print(c(B[,ac(as.numeric(t)+1)],prj(F[,t],C[,t],B[,t],r,K)))     
   B[,ac(as.numeric(t)+1)]=prj(F[,t],C[,t],B[,t],r,K)
   }
  
  return(list(F=F,B=B))}

slv<-function(C,r,K,B0){
  B=window(FLQuant(K*B0,dimnames=dimnames(C)),end=max(as.numeric(dimnames(C)$year))+1)
  F=-log(1-C/B[,-dim(B)[2]])
  for (t in rev(rev(dimnames(B)$year)[-1])){
    res=ucminf(c(F[,t]), fn = minF, #gr = gradMinF, 
               control = list(trace = -1), 
               C=C[,t],B=B[,t],r=r,K=K)
    
    F[,t]=res$par
    
    B[,ac(as.numeric(t)+1)]=prj(F[,t],C[,t],B[,t],r,K)
    }
  
  return(FLQuants(F=F,B=B))}

setGeneric('f', function(object,...) standardGeneric('f'))
setMethod( 'f',signature(object='biodyn'),
      function(object,tolVal=1e-10,niter=200){
          fCPP(catch(object),params(object)[c("r","k","b"),tolVal,niter])
        })

setGeneric('computeStock', function(object,...) standardGeneric('computeStock'))
setMethod( 'computeStock',signature(object='biodyn'),
           function(object,tolVal=1e-10,niter=200){
             stockCPP(catch(object),params(object)[c("r","k","b"),tolVal,niter])
             })


if (FALSE){
  grad(yield, .01,         B=B[,t],r=r,K=K)
  c(  gradY(  .01, C=C[,t],B=B[,t],r=r,K=K))
  
  grad(fnF,.01, C=C[,t],B=B[,t],r=r,K=K)
  c( gradF(.01, C=C[,t],B=B[,t],r=r,K=K))
  
  ggplot(mdply(data.frame(F=seq(0,.03,length.out=101)), 
               function(F) c(fnY(F,C[,t],B[,t],r,K)/dfdY(F,C[,t],B[,t],r,K))))+
    geom_line(aes(F,V1))
  
  ggplot(mdply(data.frame(F=seq(0,.7,length.out=101)), 
               function(F) c(fnY(F,C[,t],B[,t],r,K))))+
    geom_line(aes(F,V1))
  
  ggplot(mdply(data.frame(F=seq(0,.7,length.out=101)), 
               function(F) c(gradY(F,C[,t],B[,t],r,K))))+
    geom_line(aes(F,V1))
  }