R/ypr.R

Defines functions ypr

Documented in ypr

ypr<-function(age=NULL,wgt=NULL,partial=NULL,M=NULL,plus=FALSE,oldest=NULL,maxF=2,
              incrF=0.001,options=c(1),Fsol=NULL,graph=TRUE){	
  #8/21/2023 add options argument and procedures 
	if(is.null(age)) 
         stop ("age vectoris missing") 
  if(is.null(age)) 
         stop ("wgt vector is missing") 
  if(is.null(partial)) 
         stop ("partial recruitment vector is missing")
  if(is.null(M)) 
         stop ("M value or vector is missing")
  if(plus==TRUE & is.null(oldest)) stop("oldest must be specified for plus group calculation.")    
  if(any(length(age)!=c(length(age),length(wgt),length(partial))))
         stop("Length of vectors unequal")
  if(length(M)==1) M<-rep(M,length(age))
  if(length(which(!options %in% c(1,2)))>0) stop("invalid option number") 
  
  data<-as.data.frame(cbind(age,wgt,partial,M))
	YPR<-as.data.frame(cbind(rep(NA,ceiling(maxF/incrF)+1),rep(NA,ceiling(maxF/incrF)+1)))
	names(YPR)<-c("F","YPR")
  if(plus==TRUE){
       len<-oldest-min(data$age)+1
       if(oldest>max(data$age)){
          pdata<-data[rep(length(data$age),times=oldest-data$age[length(data$age)]), ] 
          pdata$age<-seq(max(data$age)+1,oldest,1)
          data<-rbind(data,pdata)
        }
    }
  if(plus==FALSE) len<-max(data$age)-min(data$age)+1
	
	if(any(options==1)){
	  F<-0.0					      
	  for(i in 1:length(YPR$F)){
            data$S<-exp(-data$partial*F-data$M) 
            data$psb[1]<-1
  		for(y in 2:len){
      	   data$psb[y]<-data$S[y-1]
    		}
         	data$psb<-cumprod(data$psb)
   		data$YPR<-((data$partial*F)/(data$partial*F+data$M))*
               (1-exp(-data$partial*F-data$M))*data$psb*data$wgt
    		YPR$YPR[i]<-sum(data$YPR)
    		YPR$F[i]<-F
    		F<-F+incrF
  	 }
   Ymax<-max(YPR$YPR)
	 Fmax<-YPR$F[which(YPR$YPR==max(YPR$YPR))]
	 s10<-((YPR$YPR[2]-YPR$YPR[1])/(YPR$F[2]-YPR$F[1]))*0.10

	F10<-Fmax/2;df<-F10/2;ok<-0;fuzz<-0.0001                
	while (ok==0){
    		  data$s1<-exp(-data$partial*F10-data$M)
    		  data$s1p<-1
    		  for(y in 2:len){
       		data$s1p[y]<-data$s1[y-1]
      }
    	data$s1p<-cumprod(data$s1p)
    	d1<-sum(((data$partial*F10)/(data$partial*F10+data$M))*
             (1-exp(-data$partial*F10-data$M))*data$s1p*data$wgt)
    	data$s2<-exp(-data$partial*(F10+0.0001)-data$M)
    	data$s2p<-1
    	for(y in 2:len){
     	    data$s2p[y]<-data$s2[y-1]
     	}
    	data$s2p<-cumprod(data$s2p)
    	d2<-sum(((data$partial*(F10+0.0001))/(data$partial*(F10+0.0001)+data$M))*
            	(1-exp(-data$partial*(F10+0.0001)-data$M))*data$s2p*data$wgt)
   	  slope<-(d2-d1)/((F10+0.0001)-F10)
   	  if(abs(s10-slope)<=fuzz) ok<-1
       if(ok==0){
        if(slope>s10) F10<-F10+df
    	  if(slope<s10) F10<-F10-df
    	  df<-df/2 
    	}
    	Y10<-d1  
	 }
}#options==1
	
#Find YPR for single F
if(any(options==2)){
    if(is.null(Fsol)) stop ("Fsol value is missing for option 2.")
      data$S<-exp(-data$partial*Fsol-data$M) 
      data$psb[1]<-1
      for(y in 2:len){
        data$psb[y]<-data$S[y-1]
      }
      data$psb<-cumprod(data$psb)
      data$YPR<-((data$partial*Fsol)/(data$partial*Fsol+data$M))*
        (1-exp(-data$partial*Fsol-data$M))*data$psb*data$wgt
    YPR2<-data.frame(F=Fsol,YPR=sum(data$YPR))
} 
	ans<-NULL
	labels<-NULL
if(any(options==1)){
  tempo<-data.frame(F=c(F10,Fmax),Yield=c(Y10,Ymax))
  rownames(tempo)<-c("F0.10","FMax")
  ans<-c(ans,list(tempo))
  labels<-c(labels,"Reference_Points")
}	
	if(any(options==2)){
	  ans<-c(ans,list(YPR2)) 
	  labels<-c(labels,"YPR_at_Fsol")
	}

	if(any(options==1)){
	  ans<-c(ans,list(YPR))
	  labels<-c(labels,"F_vs_YPR")
	}
	names(ans)<-labels  
    if(any(options==1) & graph==TRUE) {
      plot(YPR[,2]~YPR[,1],ylab="Yield-Per-Recruit",xlab="F",type="l")
    }
    return(ans)
}

Try the fishmethods package in your browser

Any scripts or data that you put into this service are public.

fishmethods documentation built on April 4, 2025, 1:47 a.m.