From 29c69309ac87ca301669ca1fceb5f6183c351cef Mon Sep 17 00:00:00 2001 From: nielshintzen Date: Tue, 5 Nov 2024 09:23:42 +0100 Subject: [PATCH] adding monteCarloStock for scanning over m --- R/monteCarloStock.R | 163 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 162 insertions(+), 1 deletion(-) diff --git a/R/monteCarloStock.R b/R/monteCarloStock.R index f3b1d2c..776cf19 100644 --- a/R/monteCarloStock.R +++ b/R/monteCarloStock.R @@ -1,4 +1,4 @@ -monteCarloStock <- function(stck,tun,sam,realisations,return.sam=FALSE,saveParsDir=NULL,set.pars=NULL,ncores=detectCores()-1, ...){ +monteCarloStock <- function(stck,tun,sam,realisations,return.sam=FALSE,saveParsDir=NULL,saveSimdataDir=NULL,set.pars=NULL,ncores=detectCores()-1, ...){ require(doParallel) ctrl <- sam@control @@ -75,6 +75,9 @@ monteCarloStock <- function(stck,tun,sam,realisations,return.sam=FALSE,saveParsD if(length(which(is.na(object$data$logobs)))>0) simdat[[i]]$logobs[which(is.na(object$data$logobs))] <- NA } + if(!is.null(saveSimdataDir)) + save(simdat,file=saveSimdataDir) + #clusterExport(cl,varlist=c("simdat","object"),envir=environment()) if(!is.null(set.pars)){ runs <- foreach(i = 1:realisations) %dopar% try(sam.fit(simdat[[i]],object$conf,object$pl,map=map,silent=T,...)) @@ -109,4 +112,162 @@ monteCarloStock <- function(stck,tun,sam,realisations,return.sam=FALSE,saveParsD ret <- list(sam=samRuns,stk=mcstck) if(!return.sam) ret <- mcstck + return(ret)} + +monteCarloStockMscan <- function(stck,tun,sam,realisations,return.sam=FALSE, + saveParsDir=NULL,saveNLLs=NULL,saveSimdataDir=NULL, + set.pars=NULL,mscan=NULL,ncores=detectCores()-1, ...){ + + require(doParallel) + ctrl <- sam@control + + #-Create new stock object with nr of realisations in iter slot + if(class(stck)=="FLStocks"){ + idxStck <- which.max(unlist(lapply(stck,function(x){x@range["maxyear"]}))) + mcstck <- propagate(stck[[idxStck]],iter=realisations) + mcstck <- window(mcstck,start=range(sam)["minyear"],end=range(sam)["maxyear"]) + mcstck@stock.n[] <- NA + mcstck@harvest[] <- NA + } else { + mcstck <- propagate(stck,iter=realisations) + mcstck <- window(mcstck,start=range(sam)["minyear"],end=range(sam)["maxyear"]) + mcstck@stock.n[] <- NA + mcstck@harvest[] <- NA + } + + #- Run first an assessment and generate new data to fit new assessment on + sam@control@simulate <- TRUE + sam@control@residuals<- FALSE + if(!is.null(set.pars)){ + if(class(stck)=="FLStock") stcks <- FLStocks(residual=stck) + dataS <- FLSAM2SAM(stcks,tun,ctrl@sumFleets,catch.vars) + confS <- ctrl2conf(ctrl,dataS) + parS <- defpar(dataS,confS) + 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.pars.internal <- merge(set.pars,matchNames,by.x="par",by.y="FLSAM") + + + mapDef <- as.list(matchNames[,"SAM"]); names(mapDef) <- matchNames[,"SAM"] + for(i in names(mapDef)) + mapDef[[i]] <- jitter(parS[[i]],factor=0.1) + if(any(duplicated(mapDef[[i]]))) + stop("mapping goes wrong, some parameters get the same starting value, retry") + map <- mapDef + + for(i in 1:nrow(set.pars.internal)){ + parS[[set.pars.internal$SAM[i]]][set.pars.internal$number[i]+1] <- set.pars.internal$value[i] + map[[set.pars.internal$SAM[[i]]]][set.pars.internal$number[i]+1] <- NA + } + map=list(logSdLogN=as.factor(map$logSdLogN), + logSdLogP=as.factor(map$logSdLogP), + logFpar=as.factor(map$logFpar), + logQpow=as.factor(map$logQpow), + logSdLogFsta=as.factor(map$logSdLogFsta), + logSdLogObs=as.factor(map$logSdLogObs)) + map <- map[which(names(map) %in% set.pars$SAM)] + + object <- FLSAM(stcks,tun,sam@control,set.pars=set.pars,return.fit=T) + } else { + object <- FLSAM(stck,tun,sam@control,return.fit=T) + } + simdat <- replicate(realisations, + c(object$data[names(object$data)!="logobs"],#all the old data + object$obj$simulate()["logobs"]),#simulated observations + simplify=FALSE) + #- Make cluster to speed up process + ncores <- ifelse(realisations0) + simdat[[i]]$logobs[which(is.na(object$data$logobs))] <- NA + } + if(!is.null(saveSimdataDir)) + save(simdat,file=saveSimdataDir) + + #Store nlls + simdatM <- simdat + storeNLL<- matrix(NA,nrow=realisations,ncol=length(mscan),dimnames=list(iter=1:realisations,addM=mscan)) + # Change input data + for(iM in mscan){ + print(iM) + for(i in 1:realisations){ + simdatM[[i]]$natMor <- simdat[[i]]$natMor + iM + } + + cl <- makeCluster(ncores) + clusterEvalQ(cl,library(stockassessment)) + registerDoParallel(cl) + + #run assessments + if(!is.null(set.pars)){ + capture.output(runs <- foreach(i = 1:realisations) %dopar% try(sam.fit(simdatM[[i]],object$conf,object$pl[-lengt(object$pl)],map=map,silent=T,newtonsteps=0)) + } else { + capture.output(runs <- foreach(i = 1:realisations) %dopar% try(sam.fit(simdatM[[i]],object$conf,object$pl[-length(object$pl)],silent=T,newtonsteps=0))) + } + + storeNLL[,ac(iM)] <- unlist(lapply(runs,function(x){try(x$opt$objective)})) + save(storeNLL,file=saveNLLs) + stopCluster(cl) + } + + #- Get the addm values, and set to default of 0 when non-convergence + addms <- mscan[apply(storeNLL,1,which.min)] + if(any(is.na(addms) | !is.finite(addms))){ + idx <- which(is.na(addms) | !is.finite(addms)) + addms[idx] <- 0 + } + + #- Final run, adjust the m vectors + for(i in 1:realisations){ + simdatM[[i]]$natMor <- simdat[[i]]$natMor + addms[i] + } + #run assessments + cl <- makeCluster(ncores) + clusterEvalQ(cl,library(stockassessment)) + registerDoParallel(cl) + if(!is.null(set.pars)){ + runs <- foreach(i = 1:realisations) %dopar% try(sam.fit(simdatM[[i]],object$conf,object$pl[-length(object$pl)],map=map,silent=T,newtonsteps=0,...)) + } else { + runs <- foreach(i = 1:realisations) %dopar% try(sam.fit(simdatM[[i]],object$conf,object$pl[-length(object$pl)],silent=T,newtonsteps=0)) + } + stopCluster(cl) #shut it down + + + #- Fill the results of the simulations + samRuns <- list() + for(i in 1:realisations){ + if(!is.na(unlist(runs[[i]]$sdrep)[1])){ + samRuns[[i]] <- SAM2FLR(runs[[i]],sam@control) + mcstck@stock.n[,,,,,i] <- samRuns[[i]]@stock.n + mcstck@harvest[,,,,,i] <- samRuns[[i]]@harvest + #mcstck@catch[,,,,,i] <- subset(params(samRuns[[i]]),name=="logCatch")$value + } + else { + samRuns[[i]] <- NA + } + } + samRuns <- as(samRuns, "FLSAMs") + + #- Save randomly drawn parameters + pars <- lapply( samRuns , function(x) params(x)$value) + random.param <- matrix(unlist(pars) , byrow= T, nrow = realisations ,dimnames =list(NULL, params(samRuns[[1]])$name )) + if(!is.null(saveParsDir)) + save(random.param,file=saveParsDir) + + if(return.sam) + ret <- list(sam=samRuns,stk=mcstck,addm=addms) + if(!return.sam) + ret <- mcstck return(ret)} \ No newline at end of file