diff --git a/PARALLEL_PROCESSING/CatchMSY_Dec2014_test.R b/PARALLEL_PROCESSING/CatchMSY_Dec2014_test.R deleted file mode 100644 index 3db4ea9..0000000 --- a/PARALLEL_PROCESSING/CatchMSY_Dec2014_test.R +++ /dev/null @@ -1,178 +0,0 @@ -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") -} diff --git a/PARALLEL_PROCESSING/CatchMSY_Dec2014_test2.R b/PARALLEL_PROCESSING/CatchMSY_Dec2014_test2.R deleted file mode 100644 index 4c72912..0000000 --- a/PARALLEL_PROCESSING/CatchMSY_Dec2014_test2.R +++ /dev/null @@ -1,178 +0,0 @@ -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_Output1.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") -} diff --git a/PARALLEL_PROCESSING/script_FAOMSY.sh b/PARALLEL_PROCESSING/script_FAOMSY.sh deleted file mode 100644 index 27ec7dd..0000000 --- a/PARALLEL_PROCESSING/script_FAOMSY.sh +++ /dev/null @@ -1,5 +0,0 @@ -#!/bin/sh -# FAOMSY -cd $1 - -java -Xmx1024M -classpath ./:./c3p0-0.9.1.2.jar:./common-configuration-scanner-1.0.1-SNAPSHOT.jar:./common-encryption-1.0.1-3.5.0.jar:./common-gcore-resources-1.2.0-3.5.0.jar:./common-gcore-stubs-1.2.0-3.5.0.jar:./common-scope-1.2.1-SNAPSHOT.jar:./common-scope-maps-1.0.2-3.5.0.jar:./commons-collections-3.1.jar:./commons-io-1.2.jar:./discovery-client-1.0.1-3.5.0.jar:./dom4j-1.6.1.jar:./ecological-engine-1.8.1-SNAPSHOT.jar:./EcologicalEngineExecutor-1.6.4-SNAPSHOT.jar:./hibernate3.jar:./ic-client-1.0.1-3.5.0.jar:./jaxen-1.1.2.jar:./jta-1.1.jar:./log4j-1.2.16.jar:./mongo-java-driver-2.12.4.jar:./postgresql-8.4-702.jdbc4.jar:./slf4j-api-1.6.0.jar:./slf4j-log4j12-1.6.0.jar:./storage-manager-core-2.1.3-3.6.0.jar:./storage-manager-wrapper-2.1.0-3.5.0.jar:./xalan-2.6.0.jar:./xpp3_min-1.1.4c.jar:./xstream-1.3.1.jar:./YASMEEN-matcher-1.2.0.1.jar:./YASMEEN-parser-1.2.0.jar org.gcube.dataanalysis.executor.nodes.algorithms.FAOMSY $2 execution.output