git-svn-id: https://svn.d4science.research-infrastructures.eu/gcube/trunk/data-analysis/EcologicalEngineSmartExecutor@114113 82a268e6-3cf1-43bd-a215-b396298e98cf
This commit is contained in:
parent
c8e72ae54b
commit
67d403919f
|
|
@ -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")
|
|
||||||
}
|
|
||||||
|
|
@ -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")
|
|
||||||
}
|
|
||||||
|
|
@ -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
|
|
||||||
Loading…
Reference in New Issue