Added combine(FLBRP) to take care of params and refpts

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: FLBRP
 Title: Reference Points for Fisheries Management
-Version: 2.5.9.9008
+Version: 2.5.9.9011
 Authors@R: c(
   person("Iago", "Mosqueira", email = "iago.mosqueira@wur.nl", role = "cre"),
   person("Laurence T.", "Kell", email = "laurie@seaplusplus.co.uk", role = "aut"),
@@ -34,6 +34,7 @@ Collate:
   fwdWindow.R
   plots.R
   spr.R
+  statistics.R
 Additional_repositories: http://flr-project.org/R
 VignetteBuilder: knitr
 LazyLoad: Yes

NAMESPACE

@@ -48,6 +48,7 @@ importFrom("utils",
   "unstack")
 
 export(
+  "abiF",
   "procerr",
   "remap")
 
@@ -64,6 +65,7 @@ exportMethods(
   "bycatch.harvest<-",
   "bycatch.wt",
   "bycatch.wt<-",
+  "combine",
   "computeRefpts",
   "discards.obs",
   "discards.obs<-",

R/coerce.R

@@ -66,7 +66,8 @@ setAs('FLBRP', 'FLStock',
     })
 
     # CREATE FLStock
-    res <- do.call("FLStock", c(full, shor, list(name=name(from), desc=desc(from))))
+    res <- do.call("FLStock", c(full, shor,
+      list(name=name(from), desc=desc(from))))
 
     # COMPUTE totals
     landings(res) <- computeLandings(res)

R/methods.R

@@ -74,9 +74,9 @@ setMethod('brp', signature(object='FLBRP'),
       refpts@.Data["virgin", "harvest",] <- 0
     }
 
-    # RUN over 500 iter blocks
+    # RUN over 250 iter blocks
     bls <- split(seq(iter), ceiling(seq_along(seq(iter)) / 250))
- 
+
     res <- lapply(bls, function(i) {
       srpars <- iter(params(object), i)
       .Call("brp", iter(object, i), 
@@ -85,10 +85,10 @@ setMethod('brp', signature(object='FLBRP'),
         PACKAGE = "FLBRP")
       }
     )
+
+    out <- Reduce(combine, res)
 
-    res <- Reduce(combine, res)
-
-    return(res)
+    return(out)
 
     res <- .Call("brp", object, refpts, SRNameCode(SRModelName(object@model)),
       FLQuant(c(params(object)),dimnames=dimnames(params(object))),
@@ -684,3 +684,18 @@ setMethod("leslie", signature(object="FLBRP",  fec="missing"),
   }
 )
 # }}}
+
+# combine {{{
+setMethod('combine', signature(x='FLBRP', y='FLBRP'),
+  function(x, y, ..., check=FALSE) {
+
+    res <- callNextMethod()
+
+    params(res) <- Reduce(combine, lapply(c(list(x,y), list(...)), params))
+
+    refpts(res) <- Reduce(combine, lapply(c(list(x,y), list(...)), refpts))
+
+    return(res)
+ }
+)
+# }}}

R/statistics.R

@@ -0,0 +1,77 @@
+# statistics.R - DESC
+# /home/mosqu003/FLR/pkgs/mine/FLBRP/R/statistics.R
+
+# Copyright (c) WUR, 2023.
+# Author: Iago MOSQUEIRA (WMR) <iago.mosqueira@wur.nl>
+#
+# Distributed under the terms of the EUPL-1.2
+
+
+# abiF {{{
+
+#' abiF
+#' Computes ABI for target F, e.g. ABImsy (Griffith et al. 2023)
+#'
+#' @param stock object of class FLStock 
+#' @param target target F at equilibrium, e.g. Fmsy
+#' @param threshold quantile ageref treshold, default 0.9
+#' @param brp An optional input FLBRP object
+#' @return FLQuant 
+#' @export
+#' @examples
+#' data(ple4)
+#' abimsy <- abiF(ple4, target=0.22, threshold=0.9)
+#' plot(abimsy) + ylim(0, 2) +
+#'  geom_hline(yintercept = 1) + ylab(expression(ABI[MSY]))
+
+abiF <- function(stock, target=0.2, threshold=0.9, brp=NULL) {
+
+  # COMPUTE FLBRP
+  if(is.null(brp)) {
+    equ <- FLBRP(stock)
+  } else if(is(brp, 'FLBRP')) {
+    equ <- brp
+    if(missing(target)) {
+      target <- fmsy(brp)
+    }
+  } else {
+    stop("brp must be of class 'FLBRP'")
+  }
+
+  # SET single F target
+  fbar(equ) <- FLQuant(target, dim=c(rep(1,5), dim(stock)[6]))
+
+  # COMPUTE for F target
+  equ <- brp(equ)
+
+  # CONVERT to FLStock
+  eqs <- as(equ, "FLStock")
+
+  # REMOVE first age
+  na <- stock.n(eqs)[-1, ]
+
+  ages <- seq(dims(na)$min, dims(na)$max)
+
+  cums <- apply(na, 2:6, cumsum)
+
+  nthresh <- quantSums(na * threshold)
+
+  maxage <- dims(stock)$max
+
+  aref <- pmin(c(apply((nthresh %-% cums) ^ 2, 2:6,
+    function(x) which(x == min(x)))),
+    range(equ)["plusgroup"] - 1)
+
+  # PROPORTION
+  rp <- unlist(Map(function(x,y) sum(x[ac(seq(y + 1, max(ages)))]) / sum(x),
+    x=split(na), y=aref))
+
+  # CALCULATE
+  res <- Reduce(combine, Map(function(x, y, z)
+    quantSums(x[ac(seq(y + 1, maxage)), ]) / quantSums(x[-1, ]) / z,
+    x=split(stock.n(stock)), y=aref, z=rp))
+
+  return(res)
+}
+
+# }}}