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

@@ -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")
+}