From dd1e570c5c1b5f71a3b2e16f0371de05c6ed23ff Mon Sep 17 00:00:00 2001 From: Gianpaolo Coro Date: Thu, 9 Apr 2015 11:07:29 +0000 Subject: [PATCH] git-svn-id: https://svn.d4science.research-infrastructures.eu/gcube/trunk/data-analysis/EcologicalEngineSmartExecutor@113924 82a268e6-3cf1-43bd-a215-b396298e98cf --- PARALLEL_PROCESSING/CatchMSY_Dec2014_test.R | 178 ++++++++++++++++++++ 1 file changed, 178 insertions(+) create mode 100644 PARALLEL_PROCESSING/CatchMSY_Dec2014_test.R diff --git a/PARALLEL_PROCESSING/CatchMSY_Dec2014_test.R b/PARALLEL_PROCESSING/CatchMSY_Dec2014_test.R new file mode 100644 index 0000000..3db4ea9 --- /dev/null +++ b/PARALLEL_PROCESSING/CatchMSY_Dec2014_test.R @@ -0,0 +1,178 @@ +set.seed(999) ## for same random sequence +#require(hacks) + +#setwd("C:/Users/Ye/Documents/Data poor fisheries/Martell Froese Method/") + +## Read Data for stock, year=yr, catch=ct, and resilience=res. Expects space delimited file with header yr ct and years in integer and catch in real with decimal point +## For example +## stock res yr ct +## cap-icel Medium 1984 1234.32 + +## filename <- "RAM_MSY.csv" +##filename <- "ICESct2.csv" + +cat("Step 1","\n") +TestRUN <- F # if it is true, just run on the test samples, false will go for a formal run! + +filename <- "D20.csv" +outfile <- "CatchMSY_Output.csv" +outfile2 <- paste("NonProcessedSpecies.csv",sep="") + +#cdat <- read.csv2(filename, header=T, dec=".") +cdat1 <- read.csv(filename) +cat("\n", "File", filename, "read successfully","\n") + +cat("Step 2","\n") +if(file.exists("cdat.RData")) +{load("cdat.RData")} else +{ + + dim(cdat1) + yrs=1950:2012 + + # to set NA as 0 + cdat1[is.na(cdat1)] <- 0 + nrow <- length(cdat1[,1]) + ndatColn <- length(cdat1[1,c(-1:-12)]) + rownames(cdat1) <- NULL + + cdat <- NULL + for(i in 1:nrow) + {#i=1 + #a <- ctotal3[i,-1] + tmp=data.frame(stock=rep(as.character(cdat1[i,"Stock_ID"]),ndatColn), + species=rep(as.character(cdat1[i,"Scientific_name"]),ndatColn), + yr=yrs,ct=unlist(c(cdat1[i,c(-1:-12)])), + res=rep(cdat1[i,"ResilienceIndex"],ndatColn)) + + cdat <- rbind(cdat,tmp) + #edit(cdat) + } +} + +StockList=unique(as.character(cdat$stock)) + +colnames(cdat) + + +#stock_id <- unique(as.character(cdat$stock)) +#?? +# stock_id <- "cod-2224" ## for selecting individual stocks +# stock=stock_id +#?? + +cat("Step 3","\n") + +## FUNCTIONS are going to be used subsequently +.schaefer <- function(theta) +{ + with(as.list(theta), { ## for all combinations of ri & ki + bt=vector() + ell = 0 ## initialize ell + J=0 #Ye + for (j in startbt) + { + if(ell == 0) + { + bt[1]=j*k*exp(rnorm(1,0, sigR)) ## set biomass in first year + for(i in 1:nyr) ## for all years in the time series + { + xt=rnorm(1,0, sigR) + bt[i+1]=(bt[i]+r*bt[i]*(1-bt[i]/k)-ct[i])*exp(xt) + ## calculate biomass as function of previous year's biomass plus net production minus catch + } + + #Bernoulli likelihood, assign 0 or 1 to each combination of r and k + ell = 0 + if(bt[nyr+1]/k>=lam1 && bt[nyr+1]/k <=lam2 && min(bt) > 0 && max(bt) <=k && bt[which(yr==interyr)]/k>=interbio[1] && bt[which(yr==interyr)]/k<=interbio[2]) + ell = 1 + J=j # Ye + } + } + return(list(ell=ell,J=J)) # Ye adding J=J + + + }) +} + +sraMSY <-function(theta, N) +{ + #This function conducts the stock reduction + #analysis for N trials + #args: + # theta - a list object containing: + # r (lower and upper bounds for r) + # k (lower and upper bounds for k) + # lambda (limits for current depletion) + + + with(as.list(theta), +{ + ri = exp(runif(N, log(r[1]), log(r[2]))) ## get N values between r[1] and r[2], assign to ri + ki = exp(runif(N, log(k[1]), log(k[2]))) ## get N values between k[1] and k[2], assing to ki + itheta=cbind(r=ri,k=ki, lam1=lambda[1],lam2=lambda[2], sigR=sigR) + ## assign ri, ki, and final biomass range to itheta + M = apply(itheta,1,.schaefer) ## call Schaefer function with parameters in itheta + i=1:N + ## prototype objective function + get.ell=function(i) M[[i]]$ell + ell = sapply(i, get.ell) + get.J=function(i) M[[i]]$J # Ye + J=sapply(i,get.J) # Ye + return(list(r=ri,k=ki, ell=ell, J=J)) # Ye adding J=J +}) +} + + +getBiomass <- function(r, k, j) +{ + BT <- NULL + bt=vector() + for (v in 1:length(r)) + { + bt[1]=j[v]*k[v]*exp(rnorm(1,0, sigR)) ## set biomass in first year + for(i in 1:nyr) ## for all years in the time series + { + xt=rnorm(1,0, sigR) + bt[i+1]=(bt[i]+r[v]*bt[i]*(1-bt[i]/k[v])-ct[i])*exp(xt) + ## calculate biomass as function of previous year's biomass plus net production minus catch + } + BT=rbind(BT, t(t(bt))) + } + return(BT) +} + +## The End of Functions section + +cat("Step 4","\n") +stockLoop <- StockList +# randomly select stocks from randomly selected 5 area codes first +if(TestRUN) +{ + set.seed(999) + AreaCodeList <- unique(cdat1$AREA_Code) + sampledAC <- sample(AreaCodeList,size=5,replace=F) + stockLoop <- cdat1[cdat1$AREA_Code %in% sampledAC,c("Stock_ID")] +} + +#setup counters +counter1 <- 0 +counter2 <- 0 + +cat("Step 4","\n") +## Loop through stocks +#write.table("x",file=outfile,append = FALSE, row.names = FALSE,col.names=FALSE,sep=",") +write.table("x",file=outfile2,append = FALSE, row.names = FALSE,col.names=FALSE,sep=",") + +for(stock in stockLoop) +{ + t0<-Sys.time() + xr <- runif(1, 1.0, 10000) + x1<-c(paste("processed",xr,sep=",")) + xr <- runif(1, 1.0, 10000) + x2<-c(paste("non processed",xr,sep=",")) + #write.table(x1,file=outfile,append = T, row.names = FALSE,col.names=FALSE,sep=",") + write.table(x2,file=outfile2,append = T, row.names = FALSE,col.names=FALSE,sep=",") + + cat("Elapsed: ",Sys.time()-t0," \n") +}