adding functionality for set.pars to simulate function

flr · Apr 25, 2024 · 5370460 · 5370460
1 parent a422eb2
commit 5370460
Show file tree

Hide file tree

Showing 2 changed files with 26 additions and 6 deletions.
diff --git a/R/SAM2FLR.R b/R/SAM2FLR.R
@@ -151,7 +151,7 @@ SAM2FLR <- function(fit,ctrl){
   res@range   <- ctrl@range
 return(res)}
 
-SIM2FLR <- function(sim,fit,ctrl){
+SIM2FLR <- function(sim,fit,ctrl,set.pars=NULL){
   resList       <-list()
   for(i in 1:nrow(sim)){
     #- Create new FLSAM object and fill with ctrl elements
@@ -167,6 +167,22 @@ SIM2FLR <- function(sim,fit,ctrl){
     #- get parameter names
     pars2match <- unique(dimnames(sim)[[2]])
     keepPars    <- names(fit$pl)[which(names(fit$pl) %in% pars2match)]
+
+    #- Fix for set.pars
+    if(!is.null(set.pars)){
+      matchNames <- matrix(c("logN.vars","logSdLogN",
+                             "logP.vars","logSdLogP",
+                             "catchabilities","logFpar",
+                             "power.law.exps","logQpow",
+                             "f.vars","logSdLogFsta",
+                             "obs.vars","logSdLogObs"),ncol=2,byrow=T,dimnames=list(1:6,c("FLSAM","SAM")))
+      if(any(!set.pars$par %in% matchNames[,"FLSAM"]))
+          warning(paste(set.pars$par[which(!set.pars$par %in% matchNames[,"FLSAM"])],"cannot be set"))
+      set.parsSAM <- merge(set.pars,matchNames,by.x="par",by.y="FLSAM")
+      for(j in 1:nrow(set.pars))
+        fit$pl[[set.parsSAM$SAM[j]]] <- fit$pl[[set.parsSAM$SAM[j]]][-(set.parsSAM$number[j]+1)]
+    }
+
     pars       <- ac(rownames(as.data.frame(unlist(fit$pl))))
     idx       <- character()
     for(iK in keepPars){

diff --git a/R/methods.R b/R/methods.R
@@ -382,20 +382,24 @@ setMethod("coefficients",signature(object="FLSAMs"),
         }
 )
 
-setGeneric("simulate", function(x,y,z,n) standardGeneric("simulate"))
+setGeneric("simulate", function(x,y,z,n,set.pars) standardGeneric("simulate"))
 
 setMethod("simulate",signature(x="FLStock",y="FLIndices",z="FLSAM.control",
-  n='numeric'),
-          function(x,y,z,n=100){
-            fit   <- FLSAM(x,y,z,return.fit=T)
+  n='numeric',set.pars="data.frame"),
+          function(x,y,z,n=100,set.pars=NULL){
+            if(!is.null(set.pars))
+              fit   <- FLSAM(x,y,z,return.fit=T,set.pars=set.pars)
+            if(is.null(set.pars))
+              fit   <- FLSAM(x,y,z,return.fit=T)
             sdrep <- sdreport(fit$obj,getJointPrecision=T)
             sigma <- as.matrix(solve(sdrep$jointPrecision))
             mu    <- c(sdrep$par.fixed,sdrep$par.random)
             sim   <- mvrnorm(n,mu=mu,Sigma=sigma)
             colnames(sim) <- rownames(sigma)
             rownames(sim) <- 1:n
 
-            sim   <- SIM2FLR(sim,fit,z)
+            sim   <- SIM2FLR(sim,fit,z,set.pars)
+
           return(sim)})