Merge pull request #2 from calbertsen/main

flr · Nov 8, 2023 · 786ee5d · 786ee5d
2 parents a1309a2 + 9b6a27d
commit 786ee5d
Show file tree

Hide file tree

Showing 3 changed files with 51 additions and 16 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,13 +1,16 @@
 Package: FLSRTMB
 Title: FLSR in TMB
-Version: 1.1.4.9012
+Version: 1.1.4.9013
 Authors@R: 
  c(
    person("Henning", "Winker", email = "[email protected]",
     role = c("aut", "cre")),
   person("European Union", role = c("fnd", "cph")),
   person("Iago", "Mosqueira", email = "[email protected]",
-    role = c("aut"), comment = c(ORCID = "0000-0002-3252-0591"))
+    role = c("aut"), comment = c(ORCID = "0000-0002-3252-0591")),
+    person(c("Christoffer","Moesgaard"), "Albertsen",
+	   email = "[email protected]", role = c("ctb"),
+	   comment = c(ORCID = "0000-0003-0088-4363"))
  )
 Description: Estimates FLR spawner recruitment relationships in TMB
 X-schema.org-keywords: stock-recruit, fisheries, flr, R, TMB, AD

diff --git a/R/FLSR.R b/R/FLSR.R
@@ -68,12 +68,15 @@ setMethod("srrTMB", signature(object="FLSR"),
   s=NULL, s.est=TRUE, s.logitsd=20, r0.pr="missing",
   lplim=0.001, uplim=0.3, plim=lplim, pmax=uplim,
   nyears=NULL, report.sR0=FALSE, inits=NULL,
-  lower=NULL, upper=NULL, SDreport=TRUE,verbose=FALSE) {
+  lower=NULL, upper=NULL, SDreport=TRUE,verbose=FALSE,
+  d.type = c("None","A")) {
 
   silent = ifelse(verbose,1,0)
 
   s.inp = s
-
+
+  d.type = match.arg(d.type)
+
   if(is.null(nyears)) nyears = dim(ssb(object))[2]
   if(is.null(plim)){
     if(is.null(s)& s.est){ plim=0.01} else {plim=0.1}
@@ -199,9 +202,9 @@ setMethod("srrTMB", signature(object="FLSR"),
 
 
     if(is.null(lower))
-      lower <- c(min(log(rec)), log(0.05),-100)
+      lower <- c(min(log(rec)), log(0.05),-100,-10)
     if(is.null(upper))
-      upper <- c(max(log(rec * 20)), log(1.5),100)
+      upper <- c(max(log(rec * 20)), log(1.5),100,10)
 
     # SET TMB input
     inp <- list(
@@ -210,17 +213,20 @@ setMethod("srrTMB", signature(object="FLSR"),
                   spr0y = spr0.yr,spr0=spr0ref,plim=plim, nyears=length(ssb),slim=ll,smax=ul,
 
                   # model
-                  Rmodel = which(model==c("bevholtSV","rickerSV","segreg"))-1),
+                  Rmodel = which(model==c("bevholtSV","rickerSV","segreg"))-1,
+                  depensationModel = which(d.type == c("None","A")) - 1),
       # inits
-      Params = list(log_r0 = inits[1], log_sigR = inits[2],logit_s=inits[3]),
+      Params = list(log_r0 = inits[1], log_sigR = inits[2],logit_s=inits[3], log_d = 0),
       # bounds
       lower=lower, upper=upper,
       #
       ReportSD = SDreport
     )
 
     # Compile TMB inputs 
-    Map = list()
+      Map = list()
+      if(d.type == "None")
+          Map$log_d = factor(NA)
     # Turn off steepness estimation
     if(!s.est) Map[["logit_s"]] = factor( NA ) 
 
@@ -240,7 +246,8 @@ setMethod("srrTMB", signature(object="FLSR"),
     Report <- Obj$report()
 
     if(SDreport) {
-      SD <- try(TMB::sdreport(Obj))
+        SD <- try(TMB::sdreport(Obj))
+        attr(object,"SD") <- SD
     }
 
     # LOAD output in FLSR
@@ -268,9 +275,9 @@ setMethod("srrTMB", signature(object="FLSR"),
 
 
     if(model!="segreg"){
-      attr(object,"SV") = data.frame(s=Report$s,sigmaR=Report$sigR,R0=Report$r0,rho=rho,B0=Report$r0*spr0ref)
+      attr(object,"SV") = data.frame(s=Report$s,sigmaR=Report$sigR,R0=Report$r0,rho=rho,B0=Report$r0*spr0ref, d = ifelse(d.type=="None",NA,Report$d))
     } else{
-      attr(object,"SV") = data.frame(s=NA,sigmaR=Report$sigR,R0=Report$r0,rho=rho,B0=Report$r0*spr0ref)
+      attr(object,"SV") = data.frame(s=NA,sigmaR=Report$sigR,R0=Report$r0,rho=rho,B0=Report$r0*spr0ref, d = ifelse(d.type=="None",NA,Report$d))
 
     }
 
@@ -281,7 +288,7 @@ setMethod("srrTMB", signature(object="FLSR"),
   } # End TMB model
 
   attr(object,"settings") = list(s=s.inp,s.est=s.est,s.logitsd=s.logitsd,spr0=spr0,lplim=lplim,uplim=uplim,nyears=nyears,
-                                 inits=inits,lower=lower,upper=upper)
+                                 inits=inits,lower=lower,upper=upper,d.type=d.type)
 
 
 

diff --git a/src/FLSRTMB.cpp b/src/FLSRTMB.cpp
@@ -25,15 +25,18 @@ Type objective_function<Type>::operator() () {
   DATA_SCALAR( slim ); 
   DATA_SCALAR( smax ); 
   DATA_INTEGER(Rmodel); // Recruitment model
+  DATA_INTEGER(depensationModel);
 
   // Parameters
   PARAMETER(log_r0);
   PARAMETER(log_sigR);
   PARAMETER(logit_s);
+  PARAMETER(log_d);
   // Derived quantities
   Type r0 = exp(log_r0);
   Type sigR = exp(log_sigR);
   Type sinit = slim+0.001 + (smax-slim-0.001)*1/(1+exp(-logit_s));
+  Type d = exp(log_d);
 
   vector<Type> log_rec_hat( nyears );
   vector<Type> vy = spr0y * r0;
@@ -49,7 +52,15 @@ Type objective_function<Type>::operator() () {
   if(Rmodel==0){ // bevholtSV()
     s = sinit;
     for( int t=0; t< nyears; t++){
-      log_rec_hat(t) = log(4.0 * s * r0 *ssb(t) / (vy(t)*(1.0-s)+ssb(t)*(5.0*s-1.0)));//-pow(sigR,2)/2.0;
+      // log_rec_hat(t) = log(4.0 * s * r0 *ssb(t) / (vy(t)*(1.0-s)+ssb(t)*(5.0*s-1.0)));//-pow(sigR,2)/2.0;
+      // Rescale SSB by point of 50% recruitment
+      Type rmaxx = 4.0 * s * r0 / (5.0 * s - 1.0);
+      Type s50x = vy(t) * (1 - s) / (5.0 * s - 1.0);
+      Type ssbx = ssb(t) / s50x;
+      if(depensationModel == 1){
+	ssbx = pow(ssbx, d);
+      }
+      log_rec_hat(t) = log(rmaxx / (1.0 + 1.0 / ssbx));
     }
   }
 
@@ -58,16 +69,29 @@ Type objective_function<Type>::operator() () {
       s = sinit; // removed *20
       b = log(5.0*s)/(0.8*vy(t));
       a = exp(b*vy(t))/spr0;
-      log_rec_hat(t) = log(a*ssb(t)*exp(-b*ssb(t)));
+      //log_rec_hat(t) = log(a*ssb(t)*exp(-b*ssb(t)));
       //log_rec_hat(t) = log(r0 * ssb(t) / v * exp(s*(1.0-ssb(t)/v)));
+      // Rescale SSB by point of maximum recruitment (1/b);
+      Type ssbx = ssb(t) * b;
+      if(depensationModel == 1){
+	ssbx = pow(ssbx, d);
+      }
+      log_rec_hat(t) = log(a / b * ssbx * exp(-ssbx));
     }
   }
 
   if(Rmodel==2){ // segreg() aka Hockey Stick
     for( int t=0; t< nyears; t++){
       s = sinit;
       //log_rec_hat(t) = log(r0)+log(2.5*s/v*(ssb(t)+0.2*v/s-pow(pow(ssb(t)-0.2*v/s,2.0),0.5)));//-pow(sigR,2)/2.0;
-      log_rec_hat(t) = log(r0)+log(0.5/plim*s/vy(t)*(ssb(t)+plim*vy(t)/s-pow(pow(ssb(t)-plim*vy(t)/s,2.0),0.5)));
+      // log_rec_hat(t) = log(r0)+log(0.5/plim*s/vy(t)*(ssb(t)+plim*vy(t)/s-pow(pow(ssb(t)-plim*vy(t)/s,2.0),0.5)));
+      // Re-parameterize to have breakpoint at 1 by scaling with breakpoint
+      Type ssbx = ssb(t) / (plim*vy(t)/s);
+      if(depensationModel == 1){ // Type A: R(S^d)
+	ssbx = pow(ssbx, d);
+      }
+      log_rec_hat(t) = log(r0) + log(0.5 * (ssbx + 1.0 - pow(pow(ssbx - 1.0, 2.0), 0.5)));
+
     }
   }
 
@@ -113,6 +137,7 @@ Type objective_function<Type>::operator() () {
   REPORT( a );
   REPORT( b );
   REPORT( s );
+  REPORT( d ) ;
   ADREPORT(log(a));
   ADREPORT( log(b));