diff --git a/PARALLEL_PROCESSING/CMSY_22_noplot.R b/PARALLEL_PROCESSING/CMSY_22_noplot.R deleted file mode 100644 index 6541ff6..0000000 --- a/PARALLEL_PROCESSING/CMSY_22_noplot.R +++ /dev/null @@ -1,696 +0,0 @@ -##-------------------------------------------------------- -## CMSY analysis with estimation of total biomass, including Bayesian Schaefer -## written by Rainer Froese with support from Gianpaolo Coro in 2013-2014 -## This version adjusts biomass to average biomass over the year -## It also contains the FutureCrash option to improve prediction of final biomass -## Version 21 adds the purple point to indicate the 25th percentile of final biomass -## Version 22 accepts that no biomass or CPUE area available -##-------------------------------------------------------- -library(R2jags) # Interface with JAGS -library(coda) - -#----------------------------------------- -# Some general settings -#----------------------------------------- -# set.seed(999) # use for comparing results between runs -rm(list=ls(all=TRUE)) # clear previous variables etc -options(digits=3) # displays all numbers with three significant digits as default -graphics.off() # close graphics windows from previous sessions - -#----------------------------------------- -# General settings for the analysis -#----------------------------------------- -sigR <- 0.02 # overall process error; 0.05 works reasonable for simulations, 0.02 for real data; 0 if deterministic model -n <- 10000 # initial number of r-k pairs -batch.mode <- T # set to TRUE to suppress graphs -write.output <- T # set to true if table of output is wanted -FutureCrash <- "No" - -#----------------------------------------- -# Start output to screen -#----------------------------------------- -cat("-------------------------------------------\n") -cat("Catch-MSY Analysis,", date(),"\n") -cat("-------------------------------------------\n") - -#------------------------------------------ -# Read data and assign to vectors -#------------------------------------------ -# filename_1 <- "AllStocks_Catch4.csv" -# filename_2 <- "AllStocks_ID4.csv" -# filename_1 <- "SimCatch.csv" -# filename_2 <- "SimSpec.csv" -# filename_2 <- "SimSpecWrongS.csv" -# filename_2 <- "SimSpecWrongI.csv" -# filename_2 <- "SimSpecWrongF.csv" -# filename_2 <- "SimSpecWrongH.csv" -# filename_2 <- "SimSpecWrongL.csv" -# filename_1 <- "FishDataLim.csv" -# filename_2 <- "FishDataLimSpec.csv" - filename_1 <- "WKLIFE4Stocks.csv" - filename_2 <- "WKLIFE4ID.csv" - -outfile<-"outfile" -outfile.txt <- "outputfile.txt" - -cdat <- read.csv(filename_1, header=T, dec=".", stringsAsFactors = FALSE) -cinfo <- read.csv(filename_2, header=T, dec=".", stringsAsFactors = FALSE) -cat("Files", filename_1, ",", filename_2, "read successfully","\n") - -# Stocks with total biomass data and catch data from StartYear to EndYear -# stocks <- sort(as.character(cinfo$stock)) # All stocks -stocks<-"HLH_M07" - -# select one stock after the other -for(stock in stocks) { - # assign data from cinfo to vectors - res <- as.character(cinfo$Resilience[cinfo$stock==stock]) - StartYear <- as.numeric(cinfo$StartYear[cinfo$stock==stock]) - EndYear <- as.numeric(cinfo$EndYear[cinfo$stock==stock]) - r_low <- as.numeric(cinfo$r_low[cinfo$stock==stock]) - r_hi <- as.numeric(cinfo$r_hi[cinfo$stock==stock]) - stb_low <- as.numeric(cinfo$stb_low[cinfo$stock==stock]) - stb_hi <- as.numeric(cinfo$stb_hi[cinfo$stock==stock]) - intyr <- as.numeric(cinfo$intyr[cinfo$stock==stock]) - intbio_low <- as.numeric(cinfo$intbio_low[cinfo$stock==stock]) - intbio_hi <- as.numeric(cinfo$intbio_hi[cinfo$stock==stock]) - endbio_low <- as.numeric(cinfo$endbio_low[cinfo$stock==stock]) - endbio_hi <- as.numeric(cinfo$endbio_hi[cinfo$stock==stock]) - Btype <- as.character(cinfo$Btype[cinfo$stock==stock]) - FutureCrash <- as.character(cinfo$FutureCrash[cinfo$stock==stock]) - comment <- as.character(cinfo$comment[cinfo$stock==stock]) - - - # extract data on stock - yr <- as.numeric(cdat$yr[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear]) - ct <- as.numeric(cdat$ct[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])/1000 ## assumes that catch is given in tonnes, transforms to '000 tonnes - if(Btype=="observed" | Btype=="CPUE" | Btype=="simulated") { - bt <- as.numeric(cdat$TB[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])/1000 ## assumes that biomass is in tonnes, transforms to '000 tonnes - } else {bt <- NA} - nyr <- length(yr) # number of years in the time series - - -if(Btype!="observed") {bio <- bt} -# change biomass to moving average as assumed by Schaefer (but not for simulations or CPUE) -# for last year use reported bio -if(Btype=="observed") { - ma <- function(x){filter(x,rep(1/2,2),sides=2)} - bio <- ma(bt) - bio[length(bio)] <- bt[length(bt)] } - - # initialize vectors for viable r, k, bt - rv.all <- vector() - kv.all <- vector() - btv.all <- matrix(data=vector(),ncol=nyr+1) - - - - #---------------------------------------------------- - # Determine initial ranges for parameters and biomass - #---------------------------------------------------- - # initial range of r from input file - if(is.na(r_low)==F & is.na(r_hi)==F) { - start_r <- c(r_low,r_hi) - } else { - # initial range of r and CatchMult values based on resilience - if(res == "High") { - start_r <- c(0.6,1.5)} else if(res == "Medium") { - start_r <- c(0.2,0.8)} else if(res == "Low") { - start_r <- c(0.05,0.5)} else { # i.e. res== "Very low" - start_r <- c(0.015,0.1)} - } - - - # initial range of k values, assuming k will always be larger than max catch - # and max catch will never be smaller than a quarter of MSY - - start_k <- c(max(ct),16*max(ct)/start_r[1]) - - # initial biomass range from input file - if(is.na(stb_low)==F & is.na(stb_hi)==F) { - startbio <- c(stb_low,stb_hi) - } else { - # us low biomass at start as default - startbio <- c(0.1,0.5) - } - - MinYear <- yr[which.min(ct)] - MaxYear <- yr[which.max(ct)] - # use year and biomass range for intermediate biomass from input file - if(is.na(intbio_low)==F & is.na(intbio_hi)==F) { - intyr <- intyr - intbio <- c(intbio_low,intbio_hi) - # else if year of minimum catch is at least 3 years away from StartYear and EndYear of series, use min catch - } else if((MinYear - StartYear) > 3 & (EndYear - MinYear) > 3 ) { - # assume that biomass range in year before minimum catch was 0.01 - 0.4 - intyr <- MinYear-1 - intbio <- c(0.01,0.4) - # else if year of max catch is at least 3 years away from StartYear and EndYear of series, use max catch - } else if((MaxYear - StartYear) > 3 & (EndYear - MaxYear) > 3 ) { - # assume that biomass range in year before maximum catch was 0.3 - 0.9 - intyr <- MaxYear-1 - intbio <- c(0.3,0.9) - } else { - # assume uninformative range 0-1 in mid-year - intyr <- as.integer(mean(c(StartYear, EndYear))) - intbio <- c(0,1) } - # end of intbio setting - - # final biomass range from input file - if(is.na(endbio_low)==F & is.na(endbio_hi)==F) { - endbio <- c(endbio_low,endbio_hi) - } else { - # else use Catch/maxCatch to estimate final biomass - endbio <- if(ct[nyr]/max(ct) > 0.5) {c(0.4,0.8)} else {c(0.01,0.4)} - } # end of final biomass setting - - - #---------------------------------------------- - # MC with Schaefer Function filtering - #---------------------------------------------- - Schaefer <- function(ri, ki, startbio, intyr, intbio, endbio, sigR, pt) { - - # if stock is not expected to crash within 3 years if last catch continues - if(FutureCrash == "No") { - yr.s <- c(yr,EndYear+1,EndYear+2,EndYear+3) - ct.s <- c(ct,ct[yr==EndYear],ct[yr==EndYear],ct[yr==EndYear]) - nyr.s <- length(yr.s) - } else{ - yr.s <- yr - ct.s <- ct - nyr.s <- nyr - } - - # create vector for initial biomasses - startbt <-seq(from =startbio[1], to=startbio[2], by = (startbio[2]-startbio[1])/10) - # create vectors for viable r, k and bt - rv <- array(-1:-1,dim=c(length(ri)*length(startbt))) #initialize array with -1. The -1 remaining after the process will be removed - kv <- array(-1:-1,dim=c(length(ri)*length(startbt))) - btv <- matrix(data=NA, nrow = (length(ri)*length(startbt)), ncol = nyr+1) - intyr.i <- which(yr.s==intyr) # get index of intermediate year - - #loop through r-k pairs - npoints = length(ri) - nstartb = length(startbt) - - for(i in 1 : npoints) { - if (i%%1000==0) - cat(".") - - # create empty vector for annual biomasses - bt <- vector() - - # loop through range of relative start biomasses - for(j in startbt) { - # set initial biomass, including process error - bt[1]=j*ki[i]*exp(rnorm(1,0, sigR)) ## set biomass in first year - - #loop through years in catch time series - for(t in 1:nyr.s) { # for all years in the time series - xt=rnorm(1,0, sigR) # set new random process error for every year - - # calculate biomass as function of previous year's biomass plus surplus production minus catch - bt[t+1]=(bt[t]+ri[i]*bt[t]*(1-bt[t]/ki[i])-ct.s[t])*exp(xt) - - # if biomass < 0.01 k or > 1.1 k, discard r-k pair - if(bt[t+1] < 0.01*ki[i] || bt[t+1] > 1.1*ki[i]) { break } # stop looping through years, go to next upper level - - if ((t+1)==intyr.i && (bt[t+1]>(intbio[2]*ki[i]) || bt[t+1]<(intbio[1]*ki[i]))) { break } #intermediate year check - - } # end of loop of years - - # if last biomass falls without expected ranges goto next r-k pair - if(t < nyr.s || bt[yr.s==EndYear] > (endbio[2]*ki[i]) || bt[yr.s==EndYear] < (endbio[1]*ki[i])) { - next } else { - # store r, k, and bt, plot point, then go to next startbt - rv[((i-1)*nstartb)+j] <- ri[i] - kv[((i-1)*nstartb)+j] <- ki[i] - btv[((i-1)*nstartb)+j,] <- bt[1:(nyr+1)]/ki[i] #substitute a row into the matrix, exclude FutureCrash years - if(pt==T) {points(x=ri[i], y=ki[i], pch=".", cex=2, col="black") - next } - } - } # end of loop of initial biomasses - } # end of loop of r-k pairs - - rv=rv[rv!=-1] - kv=kv[kv!=-1] - btv=na.omit(btv) #delete first line - - cat("\n") - return(list(rv, kv,btv)) - } # end of Schaefer function - - #------------------------------------------------------------------ - # Uniform sampling of the r-k space - #------------------------------------------------------------------ - # get random set of r and k from log space distribution - ri1 = exp(runif(n, log(start_r[1]), log(start_r[2]))) - ki1 = exp(runif(n, log(start_k[1]), log(start_k[2]))) - - #----------------------------------------------------------------- - # Plot data and progress - #----------------------------------------------------------------- - #windows(14,9) - par(mfcol=c(2,3)) - # plot catch - plot(x=yr, y=ct, ylim=c(0,1.2*max(ct)), type ="l", bty="l", main=paste(stock,"catch"), xlab="Year", - ylab="Catch", lwd=2) - points(x=yr[which.max(ct)], y=max(ct), col="red", lwd=2) - points(x=yr[which.min(ct)], y=min(ct), col="red", lwd=2) - - # plot r-k graph - plot(ri1, ki1, xlim = start_r, ylim = start_k, log="xy", xlab="r", ylab="k", main="Finding viable r-k", pch=".", cex=2, bty="l", col="lightgrey") - - #1 - Call MC-Schaefer function to preliminary explore the space without prior information - cat(stock, ": First Monte Carlo filtering of r-k space with ",n," points\n") - MCA <- Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T) - rv.all <- append(rv.all,MCA[[1]]) - kv.all <- append(kv.all,MCA[[2]]) - btv.all <- rbind(btv.all,MCA[[3]]) - #take viable r and k values - nviablepoints = length(rv.all) - cat("* Found ",nviablepoints," viable points from ",n," samples\n"); - - - #if few points were found then resample and shrink the k log space - if (nviablepoints<=1000){ - log.start_k.new <- log(start_k) - max_attempts = 3 - current_attempts = 1 - while (nviablepoints<=1000 && current_attempts<=max_attempts){ - if(nviablepoints > 0) { - log.start_k.new[1] <- mean(c(log.start_k.new[1], min(log(kv.all)))) - log.start_k.new[2] <- mean(c(log.start_k.new[2], max(log(kv.all)))) } - n.new=n*current_attempts #add more points - ri1 = exp(runif(n.new, log(start_r[1]), log(start_r[2]))) - ki1 = exp(runif(n.new, log.start_k.new[1], log.start_k.new[2])) - cat("Shrinking k space: repeating Monte Carlo in the interval [",exp(log.start_k.new[1]),",",exp(log.start_k.new[2]),"]\n") - cat("Attempt ",current_attempts," of ",max_attempts," with ",n.new," points","\n") - MCA <- Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T) - rv.all <- append(rv.all,MCA[[1]]) - kv.all <- append(kv.all,MCA[[2]]) - btv.all <- rbind(btv.all,MCA[[3]]) - nviablepoints = length(rv.all) #recalculate viable points - cat("* Found altogether",nviablepoints," viable points \n"); - current_attempts=current_attempts+1 #increment the number of attempts - } - } - - # If tip of viable r-k pairs is 'thin', do extra sampling there - gm.rv = exp(mean(log(rv.all))) - if(length(rv.all[rv.all > 0.9*start_r[2]]) < 10) { - l.sample.r <- (gm.rv + max(rv.all))/2 - cat("Final sampling in the tip area above r =",l.sample.r,"\n") - log.start_k.new <- c(log(0.8*min(kv.all)),log(max(kv.all[rv.all > l.sample.r]))) - ri1 = exp(runif(50000, log(l.sample.r), log(start_r[2]))) - ki1 = exp(runif(50000, log.start_k.new[1], log.start_k.new[2])) - MCA <- Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T) - rv.all <- append(rv.all,MCA[[1]]) - kv.all <- append(kv.all,MCA[[2]]) - btv.all <- rbind(btv.all,MCA[[3]]) - nviablepoints = length(rv.all) #recalculate viable points - cat("Found altogether", length(rv.all), "unique viable r-k pairs and biomass trajectories\n") - } - - - # ------------------------------------------------------------ - # Bayesian analysis of catch & biomass with Schaefer model - # ------------------------------------------------------------ - if(Btype == "observed" | Btype=="simulated") { - cat("Running Schaefer MCMC analysis....\n") - mcmc.burn <- as.integer(30000) - mcmc.chainLength <- as.integer(60000) # burn-in plus post-burn - mcmc.thin = 10 # to reduce autocorrelation - mcmc.chains = 3 # needs to be at least 2 for DIC - - # Parameters to be returned by JAGS - jags.save.params=c('r','k','sigma.b', 'alpha', 'sigma.r') # - - # JAGS model - Model = "model{ - # to avoid crash due to 0 values - eps<-0.01 - # set a quite narrow variation from the expected value - sigma.b <- 1/16 - tau.b <- pow(sigma.b,-2) - - Bm[1] <- log(alpha*k) - bio[1] ~ dlnorm(Bm[1],tau.b) - - - for (t in 2:nyr){ - bio[t] ~ dlnorm(Bm[t],tau.b) - Bm[t] <- log(max(bio[t-1] + r*bio[t-1]*(1 - (bio[t-1])/k) - ct[t-1], eps)) - } - - # priors - alpha ~ dunif(0.01,1) # needed for fit of first biomass - #inverse cubic root relationship between the range of viable r and the size of the search space - inverseRangeFactor <- 1/((start_r[2]-start_r[1])^1/3) - - # give sigma some variability in the inverse relationship - sigma.r ~ dunif(0.001*inverseRangeFactor,0.02*inverseRangeFactor) - tau.r <- pow(sigma.r,-2) - rm <- log((start_r[1]+start_r[2])/2) - r ~ dlnorm(rm,tau.r) - - # search in the k space from the center of the range. Allow high variability - km <- log((start_k[1]+start_k[2])/2) - tau.k <- pow(km,-2) - k ~ dlnorm(km,tau.k) - - #end model - }" - - # Write JAGS model to file - cat(Model, file="r2jags.bug") - - ### random seed - set.seed(runif(1,1,500)) # needed in JAGS - - ### run model - jags_outputs <- jags(data=c('ct','bio','nyr', 'start_r', 'start_k'), - working.directory=NULL, inits=NULL, - parameters.to.save= jags.save.params, - model.file="r2jags.bug", n.chains = mcmc.chains, - n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength, - refresh=mcmc.burn/20, ) - - # ------------------------------------------------------ - # Results from JAGS Schaefer - # ------------------------------------------------------ - r_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$r)) - k_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$k)) -## sigma_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$sigma.b)) - alpha_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$alpha)) -## sigma.r_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$sigma.r)) - - mean.log.r.jags <- mean(log(r_out)) - SD.log.r.jags <- sd(log(r_out)) - lcl.log.r.jags <- mean.log.r.jags-1.96*SD.log.r.jags - ucl.log.r.jags <- mean.log.r.jags+1.96*SD.log.r.jags - gm.r.jags <- exp(mean.log.r.jags) - lcl.r.jags <- exp(lcl.log.r.jags) - ucl.r.jags <- exp(ucl.log.r.jags) - mean.log.k.jags <- mean(log(k_out)) - SD.log.k.jags <- sd(log(k_out)) - lcl.log.k.jags <- mean.log.k.jags-1.96*SD.log.k.jags - ucl.log.k.jags <- mean.log.k.jags+1.96*SD.log.k.jags - gm.k.jags <- exp(mean.log.k.jags) - lcl.k.jags <- exp(lcl.log.k.jags) - ucl.k.jags <- exp(ucl.log.k.jags) - mean.log.MSY.jags<- mean(log(r_out)+log(k_out)-log(4)) - SD.log.MSY.jags <- sd(log(r_out)+log(k_out)-log(4)) - gm.MSY.jags <- exp(mean.log.MSY.jags) - lcl.MSY.jags <- exp(mean.log.MSY.jags-1.96*SD.log.MSY.jags) - ucl.MSY.jags <- exp(mean.log.MSY.jags+1.96*SD.log.MSY.jags) - -} # end of MCMC Schaefer loop - -#------------------------------------ -# get results from CMSY -#------------------------------------ -# get estimate of most probable r as median of mid log.r-classes above cut-off -# get remaining viable log.r and log.k -rem.log.r <- log(rv.all[rv.all > gm.rv]) -rem.log.k <- log(kv.all[rv.all>gm.rv]) -# get vectors with numbers of r and mid values in about 25 classes -hist.log.r <- hist(x=rem.log.r, breaks=25, plot=F) -log.r.counts <- hist.log.r$counts -log.r.mids <- hist.log.r$mids -# get most probable log.r as mean of mids with counts > 0 -log.r.est <- median(log.r.mids[which(log.r.counts > 0)]) -lcl.log.r <- as.numeric(quantile(x=log.r.mids[which(log.r.counts > 0)], 0.025)) -ucl.log.r <- as.numeric(quantile(x=log.r.mids[which(log.r.counts > 0)], 0.975)) -r.est <- exp(log.r.est) -lcl.r.est <- exp(lcl.log.r) -ucl.r.est <- exp(ucl.log.r) - -# do linear regression of log k ~ log r with slope fixed to -1 (from Schaefer) -reg <- lm(rem.log.k ~ 1 + offset(-1*rem.log.r)) -int.reg <- as.numeric(reg[1]) -sd.reg <- sd(resid(reg)) -se.reg <- summary(reg)$coefficients[2] -# get estimate of log(k) from y where x = log.r.est -log.k.est <- int.reg + (-1) * log.r.est -# get estimates of CL of log.k.est from y +/- SD where x = lcl.log r or ucl.log.r -lcl.log.k <- int.reg + (-1) * ucl.log.r - sd.reg -ucl.log.k <- int.reg + (-1) * lcl.log.r + sd.reg -k.est <- exp(log.k.est) -lcl.k.est <- exp(lcl.log.k) -ucl.k.est <- exp(ucl.log.k) - -# get MSY from remaining log r-k pairs -log.MSY.est <- mean(rem.log.r + rem.log.k - log(4)) -sd.log.MSY.est <- sd(rem.log.r + rem.log.k - log(4)) -lcl.log.MSY.est <- log.MSY.est - 1.96*sd.log.MSY.est -ucl.log.MSY.est <- log.MSY.est + 1.96*sd.log.MSY.est -MSY.est <- exp(log.MSY.est) -lcl.MSY.est <- exp(lcl.log.MSY.est) -ucl.MSY.est <- exp(ucl.log.MSY.est) - -# get predicted biomass vectors as median and quantiles of trajectories -median.btv <- apply(btv.all,2, median) -lastyr.bio <- median.btv[length(median.btv)-1] -nextyr.bio <- median.btv[length(median.btv)] -lcl.btv <- apply(btv.all,2, quantile, probs=0.025) -q.btv <- apply(btv.all,2, quantile, probs=0.25) -ucl.btv <- apply(btv.all,2, quantile, probs=0.975) -lcl.lastyr.bio <- lcl.btv[length(lcl.btv)-1] -ucl.lastyr.bio <- ucl.btv[length(lcl.btv)-1] -lcl.nextyr.bio <- lcl.btv[length(lcl.btv)] -ucl.nextyr.bio <- ucl.btv[length(lcl.btv)] - -# ----------------------------------------- -# Plot results -# ----------------------------------------- -# Analysis of viable r-k pairs -plot(x=rv.all, y=kv.all, xlim=start_r, - ylim=c(0.9*min(kv.all, ifelse(Btype == "observed",k_out,NA), na.rm=T), 1.1*max(kv.all)), - pch=16, col="grey",log="xy", bty="l", - xlab="r", ylab="k", main="Analysis of viable r-k") -abline(v=gm.rv, lty="dashed") - -# plot points and best estimate from full Schaefer analysis -if(Btype == "observed"|Btype=="simulated") { - # plot r-k pairs from MCMC - points(x=r_out, y=k_out, pch=16,cex=0.5) - # plot best r-k pair from MCMC - points(x=gm.r.jags, y=gm.k.jags, pch=19, col="green") - lines(x=c(lcl.r.jags, ucl.r.jags),y=c(gm.k.jags,gm.k.jags), col="green") - lines(x=c(gm.r.jags,gm.r.jags),y=c(lcl.k.jags, ucl.k.jags), col="green") -} - -# if data are from simulation, plot true r and k -if(Btype=="simulated") { - l.stock <- nchar(stock) # get length of sim stock name - r.char <- substr(stock,l.stock-1,l.stock) # get last character of sim stock name - r.sim <- NA # initialize vector for r used in simulation - if(r.char=="_H") {r.sim=1; lcl.r.sim=0.8; ucl.r.sim=1.25} else - if(r.char=="_M") {r.sim=0.5;lcl.r.sim=0.4;ucl.r.sim=0.62} else - if(r.char=="_L") {r.sim=0.25;lcl.r.sim=0.2;ucl.r.sim=0.31} else {r.sim=0.05;lcl.r.sim=0.04;ucl.r.sim=0.062} - # plot true r-k point with error bars - points(x=r.sim, y=1000, pch=19, col="red") - # add +/- 20% error bars - lines(x=c(lcl.r.sim,ucl.r.sim), y=c(1000,1000), col="red") - lines(x=c(r.sim,r.sim), y=c(800,1250), col="red") -} - -# plot blue dot for proposed r-k, with 95% CL lines -points(x=r.est, y=k.est, pch=19, col="blue") -lines(x=c(lcl.r.est, ucl.r.est),y=c(k.est,k.est), col="blue") -lines(x=c(r.est,r.est),y=c(lcl.k.est, ucl.k.est), col="blue") - -# plot biomass graph -# determine k to use for red line in b/k plot -if(Btype=="simulated") {k2use <- 1000} else - if(Btype == "observed") {k2use <- gm.k.jags} else {k2use <- k.est} -# determine hight of y-axis in plot -max.y <- max(c(bio/k2use,ucl.btv,0.6,startbio[2], intbio[2],endbio[2]),na.rm=T) - -plot(x=yr,y=median.btv[1:nyr], lwd=2, xlab="Year", ylab="Relative biomass b/k", type="l", - ylim=c(0,max.y), bty="l", main=paste("Pred. biomass vs ", Btype,sep="")) -lines(x=yr, y=lcl.btv[1:nyr],type="l") -lines(x=yr, y=ucl.btv[1:nyr],type="l") -points(x=EndYear,y=q.btv[yr==EndYear], col="purple", cex=1.5, lwd=2) -abline(h=0.5, lty="dashed") -abline(h=0.25, lty="dotted") -lines(x=c(yr[1],yr[1]), y=startbio, col="blue") -lines(x=c(intyr,intyr), y=intbio, col="blue") -lines(x=c(max(yr),max(yr)), y=endbio, col="blue") - -# if observed biomass is available, plot red biomass line -if(Btype == "observed"|Btype=="simulated") { - lines(x=yr, y=bio/k2use,type="l", col="red", lwd=1) - } - -# if CPUE data are available, scale to predicted biomass range, plot red biomass line -if(Btype == "CPUE") { - par(new=T) # prepares for new plot on top of previous - plot(x=yr, y=bio, type="l", col="red", lwd=1, - ann=F,axes=F,ylim=c(0,1.2*max(bio, na.rm=T))) # forces this plot on top of previous one - axis(4, col="red", col.axis="red") -} - -# plot yield and biomass against equilibrium surplus parabola -max.y <-max(c(ct/MSY.est,ifelse(Btype=="observed"|Btype=="simulated",ct/gm.MSY.jags,NA),1.2),na.rm=T) -# plot parabola -x=seq(from=0,to=2,by=0.001) -y=4*x-(2*x)^2 -plot(x=x, y=y, xlim=c(0,1), ylim=c(0,max.y), type="l", bty="l",xlab="Relative biomass b/k", - ylab="Catch / MSY", main="Equilibrium curve") -# plot catch against CMSY biomass estimates -points(x=median.btv[1:nyr], y=ct/MSY.est, pch=16, col="grey") -points(x=q.btv[yr==EndYear],y=ct[yr==EndYear]/MSY.est, col="purple", cex=1.5, lwd=2) -# plot catch against observed biomass or CPUE -if(Btype == "observed"|Btype=="simulated") { - points(x=bio/k2use, y=ct/gm.MSY.jags, pch=16, cex=0.5) -} - -# plot exploitation rate u against u.msy -# get u derived from predicted CMSY biomass -u.CMSY <- ct/(median.btv[1:nyr]*k.est) -u.msy.CMSY <- 1-exp(-r.est/2) # # Fmsy from CMSY expressed as exploitation rate -# get u from observed or simulated biomass -if(Btype == "observed"|Btype=="simulated") { - u.bio <- ct/bio - u.msy.bio <- 1-exp(-gm.r.jags/2) -} -# get u from CPUE -if(Btype == "CPUE") { - q=max(median.btv[1:nyr][is.na(bio)==F],na.rm=T)*k.est/max(bio,na.rm=T) - u.CPUE <- ct/(q*bio) -} - -# determine upper bound of Y-axis -max.y <- max(c(1.5, 1.2*u.CMSY/u.msy.CMSY,ct[yr==EndYear]/(q.btv[yr==EndYear]*k.est)/u.msy.CMSY, - ifelse(Btype=="observed"|Btype=="simulated",max(u.bio[is.na(u.bio)==F]/u.msy.bio),0), - na.rm=T)) -# plot u from CMSY -plot(x=yr,y=u.CMSY/u.msy.CMSY, type="l", bty="l", ylim=c(0,max.y), xlab="Year", - ylab="u / u_msy", main="Exploitation rate") -abline(h=1, lty="dashed") -points(x=EndYear,y=ct[yr==EndYear]/(q.btv[yr==EndYear]*k.est)/u.msy.CMSY, col="purple", cex=1.5, lwd=2) -# plot u from biomass -if(Btype == "observed"|Btype=="simulated") lines(x=yr, y=u.bio/u.msy.bio, col="red") -# plot u from CPUE -if(Btype == "CPUE") { - par(new=T) # prepares for new plot on top of previous - plot(x=yr, y=u.CPUE, type="l", col="red", ylim=c(0, 1.2*max(u.CPUE,na.rm=T)),ann=F,axes=F) - axis(4, col="red", col.axis="red") -} -if(batch.mode == TRUE) {dev.off()} # close plot window - -# ------------------------------------------ -# print input and results to screen -cat("---------------------------------------\n") - -cat("Species:", cinfo$ScientificName[cinfo$stock==stock], "\n") -cat("Name and region:", cinfo$EnglishName[cinfo$stock==stock], ",", cinfo$Name[cinfo$stock==stock], "\n") -cat("Stock:",stock,"\n") -cat("Catch data used from years", min(yr),"-", max(yr), "\n") -cat("Prior initial relative biomass =", startbio[1], "-", startbio[2], "\n") -cat("Prior intermediate rel. biomass=", intbio[1], "-", intbio[2], "in year", intyr, "\n") -cat("Prior final relative biomass =", endbio[1], "-", endbio[2], "\n") -cat("If current catches continue, is the stock likely to crash within 3 years?",FutureCrash,"\n") -cat("Prior range for r =", format(start_r[1],digits=2), "-", format(start_r[2],digits=2), - ", prior range for k =", start_k[1], "-", start_k[2],"\n") - -# if data are simulated, print true r-k -if(filename_1=="SimCatch.csv") { -cat("True r =", r.sim, "(because input data were simulated with Schaefer model)\n") -cat("True k = 1000 \n") -cat("True MSY =", 1000*r.sim/4,"\n") -cat("True biomass in last year =",bio[length(bio)],"or",bio[length(bio)]/1000,"k \n") -cat("True mean catch / MSY ratio =", mean(ct)/(1000*r.sim/4),"\n") -} -# print results from full Schaefer if available -if(Btype == "observed"|Btype=="simulated") { -cat("Results from Bayesian Schaefer model using catch & biomass (",Btype,")\n") -cat("MSY =", gm.MSY.jags,", 95% CL =", lcl.MSY.jags, "-", ucl.MSY.jags,"\n") -cat("Mean catch / MSY =", mean(ct)/gm.MSY.jags,"\n") -if(Btype != "CPUE") { - cat("r =", gm.r.jags,", 95% CL =", lcl.r.jags, "-", ucl.r.jags,"\n") - cat("k =", gm.k.jags,", 95% CL =", lcl.k.jags, "-", ucl.k.jags,"\n") - } -} -# results of CMSY analysis -cat("Results of CMSY analysis \n") -cat("Altogether", nviablepoints,"unique viable r-k pairs were found \n") -cat(nviablepoints-length(rem.log.r),"r-k pairs above the initial geometric mean of r =", gm.rv, "were analysed\n") -cat("r =", r.est,", 95% CL =", lcl.r.est, "-", ucl.r.est,"\n") -cat("k =", k.est,", 95% CL =", lcl.k.est, "-", ucl.k.est,"\n") -cat("MSY =", MSY.est,", 95% CL =", lcl.MSY.est, "-", ucl.MSY.est,"\n") -cat("Predicted biomass in last year =", lastyr.bio, "2.5th perc =", lcl.lastyr.bio, - "97.5th perc =", ucl.lastyr.bio,"\n") -cat("Predicted biomass in next year =", nextyr.bio, "2.5th perc =", lcl.nextyr.bio, - "97.5th perc =", ucl.nextyr.bio,"\n") -cat("----------------------------------------------------------\n") - -## Write some results into outfile -if(write.output == TRUE) { -# write data into csv file - output = data.frame(cinfo$ScientificName[cinfo$stock==stock], stock, StartYear, EndYear, mean(ct)*1000, - ifelse(Btype=="observed"|Btype=="simulate",bio[length(bio)],NA), # last biomass on record - ifelse(Btype == "observed"|Btype=="simulated",gm.MSY.jags,NA), # full Schaefer - ifelse(Btype == "observed"|Btype=="simulated",lcl.MSY.jags,NA), - ifelse(Btype == "observed"|Btype=="simulated",ucl.MSY.jags,NA), - ifelse(Btype == "observed"|Btype=="simulated",gm.r.jags,NA), - ifelse(Btype == "observed"|Btype=="simulated",lcl.r.jags,NA), - ifelse(Btype == "observed"|Btype=="simulated",ucl.r.jags,NA), - ifelse(Btype == "observed"|Btype=="simulated",gm.k.jags,NA), - ifelse(Btype == "observed"|Btype=="simulated",lcl.k.jags,NA), - ifelse(Btype == "observed"|Btype=="simulated",ucl.k.jags,NA), - r.est, lcl.r.est, ucl.r.est, # CMSY r - k.est, lcl.k.est, ucl.k.est, # CMSY k - MSY.est, lcl.MSY.est, ucl.MSY.est, # CMSY r - lastyr.bio, lcl.lastyr.bio, ucl.lastyr.bio, # last year bio - nextyr.bio, lcl.nextyr.bio, ucl.nextyr.bio)# last year + 1 bio - - write.table(output, file=outfile, append = T, sep = ",", - dec = ".", row.names = FALSE, col.names = FALSE) - -# write some text into text outfile.txt - -cat("Species:", cinfo$ScientificName[cinfo$stock==stock], "\n", - "Name:", cinfo$EnglishName[cinfo$stock==stock], "\n", - "Region:", cinfo$Name[cinfo$stock==stock], "\n", - "Stock:",stock,"\n", - "Catch data used from years", min(yr),"-", max(yr),", biomass =", Btype, "\n", - "Prior initial relative biomass =", startbio[1], "-", startbio[2], "\n", - "Prior intermediate rel. biomass=", intbio[1], "-", intbio[2], "in year", intyr, "\n", - "Prior final relative biomass =", endbio[1], "-", endbio[2], "\n", - "Future crash with current catches?", FutureCrash, "\n", - "Prior range for r =", format(start_r[1],digits=2), "-", format(start_r[2],digits=2), - ", prior range for k =", start_k[1], "-", start_k[2],"\n", - file=outfile.txt,append=T) - - if(filename_1=="SimCatch.csv") { - cat(" True r =", r.sim, "(because input data were simulated with Schaefer model)\n", - "True k = 1000, true MSY =", 1000*r.sim/4,"\n", - "True biomass in last year =",bio[length(bio)],"or",bio[length(bio)]/1000,"k \n", - "True mean catch / MSY ratio =", mean(ct)/(1000*r.sim/4),"\n", - file=outfile.txt,append=T) - } - if(Btype == "observed"|Btype=="simulated") { - cat(" Results from Bayesian Schaefer model using catch & biomass \n", - "r =", gm.r.jags,", 95% CL =", lcl.r.jags, "-", ucl.r.jags,"\n", - "k =", gm.k.jags,", 95% CL =", lcl.k.jags, "-", ucl.k.jags,"\n", - "MSY =", gm.MSY.jags,", 95% CL =", lcl.MSY.jags, "-", ucl.MSY.jags,"\n", - "Mean catch / MSY =", mean(ct)/gm.MSY.jags,"\n", - file=outfile.txt,append=T) - } - cat(" Results of CMSY analysis with altogether", nviablepoints,"unique viable r-k pairs \n", - nviablepoints-length(rem.log.r),"r-k pairs above the initial geometric mean of r =", gm.rv, "were analysed\n", - "r =", r.est,", 95% CL =", lcl.r.est, "-", ucl.r.est,"\n", - "k =", k.est,", 95% CL =", lcl.k.est, "-", ucl.k.est,"\n", - "MSY =", MSY.est,", 95% CL =", lcl.MSY.est, "-", ucl.MSY.est,"\n", - "Predicted biomass last year b/k =", lastyr.bio, "2.5th perc b/k =", lcl.lastyr.bio, - "97.5th perc b/k =", ucl.lastyr.bio,"\n", - "Precautionary 25th percentile b/k =",q.btv[yr==EndYear],"\n", - "----------------------------------------------------------\n", - file=outfile.txt,append=T) - - } - -} # end of stocks loop diff --git a/PARALLEL_PROCESSING/CatchMSY_Dec2014.R b/PARALLEL_PROCESSING/CatchMSY_Dec2014.R deleted file mode 100644 index 2c50dc7..0000000 --- a/PARALLEL_PROCESSING/CatchMSY_Dec2014.R +++ /dev/null @@ -1,435 +0,0 @@ -set.seed(999) ## for same random sequence -#require(hacks) -#13/05/2015 -#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:2013 - - # 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) -#for(i in 1:5) - -{#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) - } - save(cdat,file="cdat.RData") -} - -StockList=unique(as.character(cdat$stock)) - -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 -for(stock in stockLoop) -{ -t0<-Sys.time() -##stock = "3845" # NB only for test single loop! - ## make graph file names: - b <- with(cdat1,cdat1[Stock_ID == stock,c(1,3,5,12)]) # Stock_ID,AREA_Names,Country,"Species" - bb <- do.call(paste,b) - - yr <- cdat$yr[as.character(cdat$stock)==stock] - ct <- as.numeric(cdat$ct[as.character(cdat$stock)==stock])/1000 ## assumes that catch is given in tonnes, transforms to '000 tonnes - res <- unique(as.character(cdat$res[as.character(cdat$stock)==stock])) ## resilience from FishBase, if needed, enable in PARAMETER SECTION - nyr <- length(yr) ## number of years in the time series - - cat("\n","Stock",stock,"\n") - flush.console() - - ## PARAMETER SECTION - mvlen=3 - ma=function(x,n=mvlen){filter(x,rep(1/n,n),sides=1)} - - ## If resilience is to be used, delete ## in rows 1-4 below and set ## in row 5 below - start_r <- if(res == "Very low"){c(0.015, 0.1)}else{ - if(res == "Low") {c(0.05,0.5)}else { - if(res == "High") {c(0.6,1.5)}else {c(0.2,1)} - } - } - ## Medium, or default if no res is found - ##start_r <- c(0.5,1.5) ## disable this line if you use resilience - start_k <- c(max(ct),50*max(ct)) ## default for upper k e.g. 100 * max catch - ## startbio <- c(0.8,1) ## assumed biomass range at start of time series, as fraction of k - ##startbio <- if(ct[1]/max(ct) < 0.5) {c(0.5,0.9)} else {c(0.3,0.6)} ## use for batch processing - - ## NB: Yimin's new idea on 20Jan14 - startbio<- if(mean(ct[1:5])/max(ct) < 0.3) {c(0.6,0.95)} else { - if(mean(ct[1:5])/max(ct)>0.3&mean(ct[1:5])/max(ct)<0.6) {c(0.3,0.7)} else { - c(0.2,0.6)}} - - interyr <- yr[2] ## interim year within time series for which biomass estimate is available; set to yr[2] if no estimates are available - interbio <- c(0, 1) ## biomass range for interim year, as fraction of k; set to 0 and 1 if not available - ## finalbio <- c(0.8, 0.9) ## biomass range after last catches, as fraction of k - ## finalbio <- if(ct[nyr]/max(ct) > 0.5) {c(0.3,0.7)} else {c(0.01,0.4)} ## use for batch processing - - ## Yimin's new stuff on 10Mar14 - #######> pre-classification - - pre.clas=ct - pre.clas[pre.clas==0]=0.1 - tx=ma(as.numeric(pre.clas),n=mvlen) - Myr=which.max(tx) - Maxc=pre.clas[which.max(tx)] - - - if(Myr==1)startbio=c(0.05,0.6)else - { - if (ct[1]/Maxc>=0.5) startbio=c(0.4,0.85) - else startbio=c(0.65,0.95) - } - - if (Myr==length(yr))finalbio=c(.4,.95) else # ie from fully to overexploited - { - if (tx[length(ct)]/Maxc>=0.5) finalbio=c(.4,.85) - else finalbio=c(.05,.6) - } - - -# if (Myr==length(yr))finalbio=c(.5,.9) -# #if (Myr=0.8) finalbio=c(.4,.8) else -# # if (tx[length(ct)]/Maxc>0.5) finalbio=c(.3,.7) else finalbio=c(.05,.6)} -# # below is the last used (20 Feb) -# if (Myr0.5) finalbio=c(.2,.8) -# else finalbio=c(.05,.6) -# } - - ##############< - n <- 30000 ## number of iterations, e.g. 100000 - sigR <- 0.0 ## process error; 0 if deterministic model; 0.05 reasonable value? 0.2 is too high - - startbt <- seq(startbio[1], startbio[2], by = 0.05) ## apply range of start biomass in steps of 0.05 - parbound <- list(r = start_r, k = start_k, lambda = finalbio, sigR) - - cat("Last year =",max(yr),", last catch =",1000*ct[nyr],"\n") - cat("Resilience =",res,"\n") - cat("Process error =", sigR,"\n") - cat("Assumed initial biomass (B/k) =", startbio[1],"-", startbio[2], " k","\n") - cat("Assumed intermediate biomass (B/k) in", interyr, " =", interbio[1],"-",interbio[2]," k","\n") - cat("Assumed final biomass (B/k) =", parbound$lambda[1],"-",parbound$lambda[2]," k","\n") - cat("Initial bounds for r =", parbound$r[1], "-", parbound$r[2],"\n") - cat("Initial bounds for k =", format(1000*parbound$k[1], digits=3), "-", format(1000*parbound$k[2],digits=3),"\n") - - flush.console() - - ## MAIN - - R1 = sraMSY(parbound, n) - - ## Get statistics on r, k, MSY and determine new bounds for r and k - r1 <- R1$r[R1$ell==1] - k1 <- R1$k[R1$ell==1] - j1 <- R1$J[R1$ell==1] # Ye - msy1 <- r1*k1/4 - mean_msy1 <- exp(mean(log(msy1))) - max_k1a <- min(k1[r1<1.1*parbound$r[1]]) ## smallest k1 near initial lower bound of r - max_k1b <- max(k1[r1*k1/4=10) - { - ## set new upper bound of r to 1.2 max r1 - parbound$r[2] <- 1.2*max(r1) - ## set new lower bound for k to 0.9 min k1 and upper bound to max_k1 - parbound$k <- c(0.9 * min(k1), max_k1) - - cat("First MSY =", format(1000*mean_msy1, digits=3),"\n") - cat("First r =", format(exp(mean(log(r1))), digits=3),"\n") - cat("New upper bound for r =", format(parbound$r[2],digits=2),"\n") - cat("New range for k =", format(1000*parbound$k[1], digits=3), "-", format(1000*parbound$k[2],digits=3),"\n") - - ## Repeat analysis with new r-k bounds - R1 = sraMSY(parbound, n) - - ## Get statistics on r, k and msy - r = R1$r[R1$ell==1] - k = R1$k[R1$ell==1] - j = R1$J[R1$ell==1] # Ye - msy = r * k / 4 - mean_ln_msy = mean(log(msy)) - - ############################################################## - ##> Ye - # BT=0 - - ## - R2<-getBiomass(r, k, j) - - #R2<-R2[-1,] - runs<-rep(1:length(r), each=nyr+1) - years=rep(yr[1]:(yr[length(yr)]+1),length=length(r)*(length(yr)+1)) - - runs=t(runs) - years=t(years) - stock_id=rep(stock,length(runs)) - R3<-cbind(as.numeric(runs), as.numeric(years), stock_id, as.numeric(R2) ) - - ## changed this, as otherwise biomass is the level of the factor below - R4<-data.frame(R3, stringsAsFactors=FALSE) - names(R4)<-c("Run", "Year", "Stock","Biomass") - - Bmsy_x<-k*0.5 - Run<-c(1:length(r)) - BMSY<-cbind(Run, Bmsy_x) - R5<-merge(R4, BMSY, by="Run", all.x=T, all.y=F) - R5$B_Bmsy<-as.numeric(paste(R5$Biomass))/R5$Bmsy_x - - ### B/Bmsy calculated for all feasible combinations of r,K,B0 - R6<-aggregate(log(B_Bmsy)~as.numeric(Year)+Stock, data=R5, - FUN=function(z){c(mean=mean(z),sd=sd(z),upr=exp(quantile(z, p=0.975)), - lwr=exp(quantile(z, p=0.025)), lwrQ=exp(quantile(z, p=0.25)), - uprQ=exp(quantile(z, p=0.75)))}) # from directly calculated from R5 becasue B_Bmsy has a lognormal dist - - R6<-data.frame(cbind(R6[,1:2],R6[,3][,1],R6[,3][,2],R6[,3][,3],R6[,3][,4],R6[,3][,5], R6[,3][,6])) - names(R6)<-c("Year", "Stock", "BoverBmsy", "BoverBmsySD","BoverBmsyUpper","BoverBmsyLower","BoverBmsylwrQ","BoverBmsyuprQ") - ##remove last entry as it is 1 greater than number of years - ## removed final year here for ease of dataframe output below - R6<-R6[-length(R6),] - ## geometric mean - GM_B_Bmsy<-exp(R6$BoverBmsy) - GM_B_BmsySD=R6$BoverBmsySD #add - ## arithmetic mean - M_B_Bmsy<-exp(R6$BoverBmsy+R6$BoverBmsySD^2/2) - - ### r,k, and MSY - - #del GM_B_Bmsy=c(rep(0,(min(yr)-1940)),GM_B_Bmsy) - #del GM_B_BmsySD=c(rep(0,(min(yr)-1940)),GM_B_BmsySD) ###### - #del M_B_Bmsy=c(rep(0,(min(yr)-1940)),M_B_Bmsy) - #del yr1=seq(1940,max(yr)) - - yr1=yr #add - - stockInfo <- with(cdat1,cdat1[Stock_ID==stock,1:12]) - temp=c(startbio[1],startbio[2],finalbio[1],finalbio[2],res, - mean(log(r)),sd(log(r)),mean(log(k)),sd(log(k)),mean(log(msy)), - sd(log(msy)),sigR,min(yr),max(yr),max(ct),length(r),GM_B_Bmsy,GM_B_BmsySD,M_B_Bmsy) - - #add, adding "GM_B_BmsySD" in the line above - - output=as.data.frame(matrix(temp,nrow=1)) - output <- cbind(stockInfo,output) - names(output) <- c(names(cdat1)[1:12],"startbio[1]","startbio[2]","finalbio[1]","finalbio[2]", - "res","mean(log(r))","sd(log(r))","mean(log(k))","sd(log(k))", - "mean(log(msy))","sd(log(msy))","sigR","min(yr)","max(yr)","max(ct)", - "length(r)",paste("GM_B_msy",yr1,sep="_"),paste("GM_B_msySD",yr1,sep="_"),paste("M_B_Bmsy",yr1,sep="_")) - - #add, adding "paste("GM_B_msySD",yr1,sep="_")"in the line above - - ######< Ye - ######################################################## - - ## plot MSY over catch data - pdf(paste(bb,"graph.pdf",sep="_")) - - par(mfcol=c(2,3)) - plot(yr, ct, type="l", ylim = c(0, max(ct)), xlab = "Year", - ylab = "Catch (1000 t)",main = paste("StockID",stock,sep=":")) - abline(h=exp(mean(log(msy))),col="red", lwd=2) - abline(h=exp(mean_ln_msy - 2 * sd(log(msy))),col="red") - abline(h=exp(mean_ln_msy + 2 * sd(log(msy))),col="red") - - hist(r, freq=F, xlim=c(0, 1.2 * max(r)), main = "") - abline(v=exp(mean(log(r))),col="red",lwd=2) - abline(v=exp(mean(log(r))-2*sd(log(r))),col="red") - abline(v=exp(mean(log(r))+2*sd(log(r))),col="red") - - plot(r1, k1, xlim = start_r, ylim = start_k, xlab="r", ylab="k (1000t)") - - hist(k, freq=F, xlim=c(0, 1.2 * max(k)), xlab="k (1000t)", main = "") - abline(v=exp(mean(log(k))),col="red", lwd=2) - abline(v=exp(mean(log(k))-2*sd(log(k))),col="red") - abline(v=exp(mean(log(k))+2*sd(log(k))),col="red") - - - plot(log(r), log(k),xlab="ln(r)",ylab="ln(k)") - abline(v=mean(log(r))) - abline(h=mean(log(k))) - abline(mean(log(msy))+log(4),-1, col="red",lwd=2) - abline(mean(log(msy))-2*sd(log(msy))+log(4),-1, col="red") - abline(mean(log(msy))+2*sd(log(msy))+log(4),-1, col="red") - - hist(msy, freq=F, xlim=c(0, 1.2 * max(msy)), xlab="MSY (1000t)",main = "") - abline(v=exp(mean(log(msy))),col="red", lwd=2) - abline(v=exp(mean_ln_msy - 2 * sd(log(msy))),col="red") - abline(v=exp(mean_ln_msy + 2 * sd(log(msy))),col="red") - - graphics.off() - - cat("Possible combinations = ", length(r),"\n") - cat("geom. mean r =", format(exp(mean(log(r))),digits=3), "\n") - cat("r +/- 2 SD =", format(exp(mean(log(r))-2*sd(log(r))),digits=3),"-",format(exp(mean(log(r))+2*sd(log(r))),digits=3), "\n") - cat("geom. mean k =", format(1000*exp(mean(log(k))),digits=3), "\n") - cat("k +/- 2 SD =", format(1000*exp(mean(log(k))-2*sd(log(k))),digits=3),"-",format(1000*exp(mean(log(k))+2*sd(log(k))),digits=3), "\n") - cat("geom. mean MSY =", format(1000*exp(mean(log(msy))),digits=3),"\n") - cat("MSY +/- 2 SD =", format(1000*exp(mean_ln_msy - 2 * sd(log(msy))),digits=3), "-", format(1000*exp(mean_ln_msy + 2 * sd(log(msy))),digits=3), "\n") - - ## Write results into outfile, in append mode (no header in file, existing files will be continued) - ## output = data.frame(stock, sigR, startbio[1], startbio[2], interbio[1], interbio[2], finalbio[1], finalbio[2], min(yr), max(yr), res, max(ct), ct[1], ct[nyr], length(r), exp(mean(log(r))), sd(log(r)), min(r), quantile(r,0.05), quantile(r,0.25), median(r), quantile(r,0.75), quantile(r,0.95), max(r), exp(mean(log(k))), sd(log(k)), min(k), quantile(k, 0.05), quantile(k, 0.25), median(k), quantile(k, 0.75), quantile(k, 0.95), max(k), exp(mean(log(msy))), sd(log(msy)), min(msy), quantile(msy, 0.05), quantile(msy, 0.25), median(msy), quantile(msy, 0.75), quantile(msy, 0.95), max(msy)) - - #write.table(output, file = outfile, append = TRUE, sep = ";", dec = ".", row.names = FALSE, col.names = FALSE) - appendPar <- ifelse(counter2==0,F,T) - colnamePar <- ifelse(counter2==0,T,F) - write.table(output, file = outfile, append = appendPar, sep = ",", dec = ".", - row.names = FALSE, col.names = colnamePar) - - counter2 <- counter2 + 1 - - } -cat("Elapsed: ",Sys.time()-t0," \n") -} ## End of stock loop, get next stock or exit diff --git a/PARALLEL_PROCESSING/CatchMSY_Dec2014_1.R b/PARALLEL_PROCESSING/CatchMSY_Dec2014_1.R deleted file mode 100644 index 51e6329..0000000 --- a/PARALLEL_PROCESSING/CatchMSY_Dec2014_1.R +++ /dev/null @@ -1,440 +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 -for(stock in stockLoop) -{ - t0<-Sys.time() -##stock = "3845" # NB only for test single loop! - ## make graph file names: - b <- with(cdat1,cdat1[Stock_ID == stock,c(1,3,5,12)]) # Stock_ID,AREA_Names,Country,"Species" - bb <- do.call(paste,b) - - yr <- cdat$yr[as.character(cdat$stock)==stock] - ct <- as.numeric(cdat$ct[as.character(cdat$stock)==stock])/1000 ## assumes that catch is given in tonnes, transforms to '000 tonnes - res <- unique(as.character(cdat$res[as.character(cdat$stock)==stock])) ## resilience from FishBase, if needed, enable in PARAMETER SECTION - nyr <- length(yr) ## number of years in the time series - - cat("\n","Stock",stock,"\n") - flush.console() - - ## PARAMETER SECTION - mvlen=3 - ma=function(x,n=mvlen){filter(x,rep(1/n,n),sides=1)} - - ## If resilience is to be used, delete ## in rows 1-4 below and set ## in row 5 below - start_r <- if(res == "Very low"){c(0.015, 0.1)}else{ - if(res == "Low") {c(0.05,0.5)}else { - if(res == "High") {c(0.6,1.5)}else {c(0.2,1)} - } - } - ## Medium, or default if no res is found - ##start_r <- c(0.5,1.5) ## disable this line if you use resilience - start_k <- c(max(ct),50*max(ct)) ## default for upper k e.g. 100 * max catch - ## startbio <- c(0.8,1) ## assumed biomass range at start of time series, as fraction of k - ##startbio <- if(ct[1]/max(ct) < 0.5) {c(0.5,0.9)} else {c(0.3,0.6)} ## use for batch processing - - ## NB: Yimin's new idea on 20Jan14 - startbio<- if(mean(ct[1:5])/max(ct) < 0.3) {c(0.6,0.95)} else { - if(mean(ct[1:5])/max(ct)>0.3&mean(ct[1:5])/max(ct)<0.6) {c(0.3,0.7)} else { - c(0.2,0.6)}} - - interyr <- yr[2] ## interim year within time series for which biomass estimate is available; set to yr[2] if no estimates are available - interbio <- c(0, 1) ## biomass range for interim year, as fraction of k; set to 0 and 1 if not available - ## finalbio <- c(0.8, 0.9) ## biomass range after last catches, as fraction of k - ## finalbio <- if(ct[nyr]/max(ct) > 0.5) {c(0.3,0.7)} else {c(0.01,0.4)} ## use for batch processing - - ## Yimin's new stuff on 10Mar14 - #######> pre-classification - - pre.clas=ct - pre.clas[pre.clas==0]=0.1 - tx=ma(as.numeric(pre.clas),n=mvlen) - Myr=which.max(tx) - Maxc=pre.clas[which.max(tx)] - - - if(Myr==1)startbio=c(0.05,0.6)else - { - if (ct[1]/Maxc>=0.5) startbio=c(0.4,0.85) - else startbio=c(0.65,0.95) - } - - if (Myr==length(yr))finalbio=c(.4,.95) else # ie from fully to overexploited - { - if (tx[length(ct)]/Maxc>=0.5) finalbio=c(.4,.85) - else finalbio=c(.05,.6) - } - - -# if (Myr==length(yr))finalbio=c(.5,.9) -# #if (Myr=0.8) finalbio=c(.4,.8) else -# # if (tx[length(ct)]/Maxc>0.5) finalbio=c(.3,.7) else finalbio=c(.05,.6)} -# # below is the last used (20 Feb) -# if (Myr0.5) finalbio=c(.2,.8) -# else finalbio=c(.05,.6) -# } - - ##############< - n <- 30000 ## number of iterations, e.g. 100000 - sigR <- 0.0 ## process error; 0 if deterministic model; 0.05 reasonable value? 0.2 is too high - - startbt <- seq(startbio[1], startbio[2], by = 0.05) ## apply range of start biomass in steps of 0.05 - parbound <- list(r = start_r, k = start_k, lambda = finalbio, sigR) - - cat("Last year =",max(yr),", last catch =",1000*ct[nyr],"\n") - cat("Resilience =",res,"\n") - cat("Process error =", sigR,"\n") - cat("Assumed initial biomass (B/k) =", startbio[1],"-", startbio[2], " k","\n") - cat("Assumed intermediate biomass (B/k) in", interyr, " =", interbio[1],"-",interbio[2]," k","\n") - cat("Assumed final biomass (B/k) =", parbound$lambda[1],"-",parbound$lambda[2]," k","\n") - cat("Initial bounds for r =", parbound$r[1], "-", parbound$r[2],"\n") - cat("Initial bounds for k =", format(1000*parbound$k[1], digits=3), "-", format(1000*parbound$k[2],digits=3),"\n") - - flush.console() - - ## MAIN - - R1 = sraMSY(parbound, n) - - ## Get statistics on r, k, MSY and determine new bounds for r and k - r1 <- R1$r[R1$ell==1] - k1 <- R1$k[R1$ell==1] - j1 <- R1$J[R1$ell==1] # Ye - msy1 <- r1*k1/4 - mean_msy1 <- exp(mean(log(msy1))) - max_k1a <- min(k1[r1<1.1*parbound$r[1]]) ## smallest k1 near initial lower bound of r - max_k1b <- max(k1[r1*k1/4=10) - { - ## set new upper bound of r to 1.2 max r1 - parbound$r[2] <- 1.2*max(r1) - ## set new lower bound for k to 0.9 min k1 and upper bound to max_k1 - parbound$k <- c(0.9 * min(k1), max_k1) - - cat("First MSY =", format(1000*mean_msy1, digits=3),"\n") - cat("First r =", format(exp(mean(log(r1))), digits=3),"\n") - cat("New upper bound for r =", format(parbound$r[2],digits=2),"\n") - cat("New range for k =", format(1000*parbound$k[1], digits=3), "-", format(1000*parbound$k[2],digits=3),"\n") - - ## Repeat analysis with new r-k bounds - R1 = sraMSY(parbound, n) - - ## Get statistics on r, k and msy - r = R1$r[R1$ell==1] - k = R1$k[R1$ell==1] - j = R1$J[R1$ell==1] # Ye - msy = r * k / 4 - mean_ln_msy = mean(log(msy)) - - ############################################################## - ##> Ye - # BT=0 - - ## - R2<-getBiomass(r, k, j) - - #R2<-R2[-1,] - runs<-rep(1:length(r), each=nyr+1) - years=rep(yr[1]:(yr[length(yr)]+1),length=length(r)*(length(yr)+1)) - - runs=t(runs) - years=t(years) - stock_id=rep(stock,length(runs)) - R3<-cbind(as.numeric(runs), as.numeric(years), stock_id, as.numeric(R2) ) - - ## changed this, as otherwise biomass is the level of the factor below - R4<-data.frame(R3, stringsAsFactors=FALSE) - names(R4)<-c("Run", "Year", "Stock","Biomass") - - Bmsy_x<-k*0.5 - Run<-c(1:length(r)) - BMSY<-cbind(Run, Bmsy_x) - R5<-merge(R4, BMSY, by="Run", all.x=T, all.y=F) - R5$B_Bmsy<-as.numeric(paste(R5$Biomass))/R5$Bmsy_x - - ### B/Bmsy calculated for all feasible combinations of r,K,B0 - R6<-aggregate(log(B_Bmsy)~as.numeric(Year)+Stock, data=R5, - FUN=function(z){c(mean=mean(z),sd=sd(z),upr=exp(quantile(z, p=0.975)), - lwr=exp(quantile(z, p=0.025)), lwrQ=exp(quantile(z, p=0.25)), - uprQ=exp(quantile(z, p=0.75)))}) # from directly calculated from R5 becasue B_Bmsy has a lognormal dist - - R6<-data.frame(cbind(R6[,1:2],R6[,3][,1],R6[,3][,2],R6[,3][,3],R6[,3][,4],R6[,3][,5], R6[,3][,6])) - names(R6)<-c("Year", "Stock", "BoverBmsy", "BoverBmsySD","BoverBmsyUpper","BoverBmsyLower","BoverBmsylwrQ","BoverBmsyuprQ") - ##remove last entry as it is 1 greater than number of years - ## removed final year here for ease of dataframe output below - R6<-R6[-length(R6),] - ## geometric mean - GM_B_Bmsy<-exp(R6$BoverBmsy) - GM_B_BmsySD=R6$BoverBmsySD #add - ## arithmetic mean - M_B_Bmsy<-exp(R6$BoverBmsy+R6$BoverBmsySD^2/2) - - ### r,k, and MSY - - #del GM_B_Bmsy=c(rep(0,(min(yr)-1940)),GM_B_Bmsy) - #del GM_B_BmsySD=c(rep(0,(min(yr)-1940)),GM_B_BmsySD) ###### - #del M_B_Bmsy=c(rep(0,(min(yr)-1940)),M_B_Bmsy) - #del yr1=seq(1940,max(yr)) - - yr1=yr #add - - stockInfo <- with(cdat1,cdat1[Stock_ID==stock,1:12]) - temp=c(startbio[1],startbio[2],finalbio[1],finalbio[2],res, - mean(log(r)),sd(log(r)),mean(log(k)),sd(log(k)),mean(log(msy)), - sd(log(msy)),sigR,min(yr),max(yr),max(ct),length(r),GM_B_Bmsy,GM_B_BmsySD,M_B_Bmsy) - - #add, adding "GM_B_BmsySD" in the line above - - output=as.data.frame(matrix(temp,nrow=1)) - output <- cbind(stockInfo,output) - names(output) <- c(names(cdat1)[1:12],"startbio[1]","startbio[2]","finalbio[1]","finalbio[2]", - "res","mean(log(r))","sd(log(r))","mean(log(k))","sd(log(k))", - "mean(log(msy))","sd(log(msy))","sigR","min(yr)","max(yr)","max(ct)", - "length(r)",paste("GM_B_msy",yr1,sep="_"),paste("GM_B_msySD",yr1,sep="_"),paste("M_B_Bmsy",yr1,sep="_")) - - #add, adding "paste("GM_B_msySD",yr1,sep="_")"in the line above - - ######< Ye - ######################################################## - - ## plot MSY over catch data - pdf(paste(bb,"graph.pdf",sep="_")) - - par(mfcol=c(2,3)) - plot(yr, ct, type="l", ylim = c(0, max(ct)), xlab = "Year", - ylab = "Catch (1000 t)",main = paste("StockID",stock,sep=":")) - abline(h=exp(mean(log(msy))),col="red", lwd=2) - abline(h=exp(mean_ln_msy - 2 * sd(log(msy))),col="red") - abline(h=exp(mean_ln_msy + 2 * sd(log(msy))),col="red") - - hist(r, freq=F, xlim=c(0, 1.2 * max(r)), main = "") - abline(v=exp(mean(log(r))),col="red",lwd=2) - abline(v=exp(mean(log(r))-2*sd(log(r))),col="red") - abline(v=exp(mean(log(r))+2*sd(log(r))),col="red") - - plot(r1, k1, xlim = start_r, ylim = start_k, xlab="r", ylab="k (1000t)") - - hist(k, freq=F, xlim=c(0, 1.2 * max(k)), xlab="k (1000t)", main = "") - abline(v=exp(mean(log(k))),col="red", lwd=2) - abline(v=exp(mean(log(k))-2*sd(log(k))),col="red") - abline(v=exp(mean(log(k))+2*sd(log(k))),col="red") - - - plot(log(r), log(k),xlab="ln(r)",ylab="ln(k)") - abline(v=mean(log(r))) - abline(h=mean(log(k))) - abline(mean(log(msy))+log(4),-1, col="red",lwd=2) - abline(mean(log(msy))-2*sd(log(msy))+log(4),-1, col="red") - abline(mean(log(msy))+2*sd(log(msy))+log(4),-1, col="red") - - hist(msy, freq=F, xlim=c(0, 1.2 * max(msy)), xlab="MSY (1000t)",main = "") - abline(v=exp(mean(log(msy))),col="red", lwd=2) - abline(v=exp(mean_ln_msy - 2 * sd(log(msy))),col="red") - abline(v=exp(mean_ln_msy + 2 * sd(log(msy))),col="red") - - graphics.off() - - - cat("Possible combinations = ", length(r),"\n") - cat("geom. mean r =", format(exp(mean(log(r))),digits=3), "\n") - cat("r +/- 2 SD =", format(exp(mean(log(r))-2*sd(log(r))),digits=3),"-",format(exp(mean(log(r))+2*sd(log(r))),digits=3), "\n") - cat("geom. mean k =", format(1000*exp(mean(log(k))),digits=3), "\n") - cat("k +/- 2 SD =", format(1000*exp(mean(log(k))-2*sd(log(k))),digits=3),"-",format(1000*exp(mean(log(k))+2*sd(log(k))),digits=3), "\n") - cat("geom. mean MSY =", format(1000*exp(mean(log(msy))),digits=3),"\n") - cat("MSY +/- 2 SD =", format(1000*exp(mean_ln_msy - 2 * sd(log(msy))),digits=3), "-", format(1000*exp(mean_ln_msy + 2 * sd(log(msy))),digits=3), "\n") - - ## Write results into outfile, in append mode (no header in file, existing files will be continued) - ## output = data.frame(stock, sigR, startbio[1], startbio[2], interbio[1], interbio[2], finalbio[1], finalbio[2], min(yr), max(yr), res, max(ct), ct[1], ct[nyr], length(r), exp(mean(log(r))), sd(log(r)), min(r), quantile(r,0.05), quantile(r,0.25), median(r), quantile(r,0.75), quantile(r,0.95), max(r), exp(mean(log(k))), sd(log(k)), min(k), quantile(k, 0.05), quantile(k, 0.25), median(k), quantile(k, 0.75), quantile(k, 0.95), max(k), exp(mean(log(msy))), sd(log(msy)), min(msy), quantile(msy, 0.05), quantile(msy, 0.25), median(msy), quantile(msy, 0.75), quantile(msy, 0.95), max(msy)) - - #write.table(output, file = outfile, append = TRUE, sep = ";", dec = ".", row.names = FALSE, col.names = FALSE) - appendPar <- ifelse(counter2==0,F,T) - colnamePar <- ifelse(counter2==0,T,F) - write.table(output, file = outfile, append = appendPar, sep = ",", dec = ".", - row.names = FALSE, col.names = colnamePar) - - counter2 <- counter2 + 1 - - } -cat("Elapsed: ",Sys.time()-t0," \n") -} ## End of stock loop, get next stock or exit diff --git a/PARALLEL_PROCESSING/DestinationDBHibernate.cfg.xml b/PARALLEL_PROCESSING/DestinationDBHibernate.cfg.xml deleted file mode 100644 index 99feeae..0000000 --- a/PARALLEL_PROCESSING/DestinationDBHibernate.cfg.xml +++ /dev/null @@ -1,17 +0,0 @@ - - - - org.postgresql.Driver - org.hibernate.connection.C3P0ConnectionProvider - jdbc:postgresql://localhost/testdb - gcube - d4science2 - org.hibernate.dialect.PostgreSQLDialect - org.hibernate.transaction.JDBCTransactionFactory - 0 - 1 - 0 - 1 - thread - - \ No newline at end of file diff --git a/PARALLEL_PROCESSING/UpdateLWR_4.R b/PARALLEL_PROCESSING/UpdateLWR_4.R deleted file mode 100644 index 9eb7a97..0000000 --- a/PARALLEL_PROCESSING/UpdateLWR_4.R +++ /dev/null @@ -1,530 +0,0 @@ -#### R and JAGS code for estimating LWR-parameters from previous studies -#### Meant for updating the ESTIMATE table in FishBase -#### Created by Rainer Froese in March 2013, including JAGS models by James Thorston -#### Modified in June 2013 to include subfamilies - -rm(list=ls(all=TRUE)) # remove previous variables and data -options(digits=3) # 3 significant digits as default -library(R2jags) # Interface with JAGS -runif(1) # sets random seed - -#### Read in data -DataFile = "RF_LWR2.csv" # RF_LWR4 was extracted from FishBase in June 2013 -Data = read.csv(DataFile, header=TRUE) -cat("Start", date(), "\n") -cat("Data file =", DataFile, "\n") -# Get unique, sorted list of Families -Fam.All <- sort(unique(as.character(Data$Family))) -Families <- Fam.All[Fam.All== "Acanthuridae" | Fam.All == "Achiridae"] - -OutFile = "LWR_Test1.csv" -JAGSFILE = "dmnorm_0.bug" - -# Get unique, sorted list of body shapes -Bshape <- sort(unique(as.character(Data$BodyShapeI))) - -#------------------------------------------ -# Functions -#------------------------------------------ - -#--------------------------------------------------------- -# Function to get the priors for the respective body shape -#--------------------------------------------------------- -Get.BS.pr <- function(BS) { - ### Assignment of priors based on available body shape information - # priors derived from 5150 LWR studies in FishBase 02/2013 - - if (BS == "eel-like") { # eel-like prior for log(a) and b - prior_mean_log10a = -2.99 - prior_sd_log10a = 0.175 - prior_tau_log10a = 1/prior_sd_log10a^2 - prior_mean_b = 3.06 - prior_sd_b = 0.0896 - prior_tau_b = 1/prior_sd_b^2 - } else - if (BS == "elongated") { # elongate prior for log(a) and b - prior_mean_log10a = -2.41 - prior_sd_log10a = 0.171 - prior_tau_log10a = 1/prior_sd_log10a^2 - prior_mean_b = 3.12 - prior_sd_b = 0.09 - prior_tau_b = 1/prior_sd_b^2 - } else - if (BS == "fusiform / normal") { # fusiform / normal prior for log(a) and b - prior_mean_log10a = -1.95 - prior_sd_log10a = 0.173 - prior_tau_log10a = 1/prior_sd_log10a^2 - prior_mean_b = 3.04 - prior_sd_b = 0.0857 - prior_tau_b = 1/prior_sd_b^2 - } else - if (BS == "short and / or deep") { # short and / or deep prior for log(a) and b - prior_mean_log10a = -1.7 - prior_sd_log10a = 0.175 - prior_tau_log10a = 1/prior_sd_log10a^2 - prior_mean_b = 3.01 - prior_sd_b = 0.0905 - prior_tau_b = 1/prior_sd_b^2 - } else - # priors across all shapes, used for missing or other BS - { - prior_mean_log10a = -2.0 - prior_sd_log10a = 0.313 - prior_tau_log10a = 1/prior_sd_log10a^2 - prior_mean_b = 3.04 - prior_sd_b = 0.119 - prior_tau_b = 1/prior_sd_b^2 - } - - # Priors for measurement error (= sigma) based on 5150 studies - # given here as shape mu and rate r, for gamma distribution - SD_rObs_log10a = 6520 - SD_muObs_log10a = 25076 - SD_rObs_b = 6808 - SD_muObs_b = 37001 - # Priors for between species variability (= sigma) based on 5150 studies for 1821 species - SD_rGS_log10a = 1372 - SD_muGS_log10a = 7933 - SD_rGS_b = 572 - SD_muGS_b = 6498 - - prior.list <- list(mean_log10a=prior_mean_log10a, sd_log10a=prior_sd_log10a, - tau_log10a=prior_tau_log10a, mean_b=prior_mean_b, sd_b=prior_sd_b, - tau_b=prior_tau_b, SD_rObs_log10a=SD_rObs_log10a, SD_muObs_log10a=SD_muObs_log10a, - SD_rObs_b=SD_rObs_b, SD_muObs_b=SD_muObs_b, SD_rGS_log10a=SD_rGS_log10a, - SD_muGS_log10a=SD_muGS_log10a, SD_rGS_b=SD_rGS_b, SD_muGS_b=SD_muGS_b) -return(prior.list) -} - -#-------------------------------------------------------------------- -# Function to do a Bayesian analysis including LWR from relatives -#-------------------------------------------------------------------- -SpecRelLWR <- function(a, b, wts, GenusSpecies, Nspecies, prior_mean_b, prior_tau_b, - prior_mean_log10a, prior_tau_log10a, SD_rObs_log10a, SD_muObs_log10a, - SD_rObs_b, SD_muObs_b, SD_rGS_log10a, SD_muGS_log10a, - SD_rGS_b, SD_muGS_b){ - ### Define JAGS model - Model = " -model { - #### Process model -- effects of taxonomy - # given the likelihood distributions and the priors, - # create normal posterior distributions for log10a, b, - # and for the process error (=between species variability sigmaGS) - - abTrue[1] ~ dnorm(prior_mean_log10a,prior_tau_log10a) - abTrue[2] ~ dnorm(prior_mean_b,prior_tau_b) - sigmaGSlog10a ~ dgamma( SD_rGS_log10a, SD_muGS_log10a) - sigmaGSb ~ dgamma( SD_rGS_b, SD_muGS_b) - - # given the posterior distributions and the process errors, - # establish for every species the expected witin-species - # parameter distributions; no correlation roGS between species - - roGS <- 0 - tauGenusSpecies[1] <- pow(sigmaGSlog10a,-2) - tauGenusSpecies[2] <- pow(sigmaGSb,-2) - for(k in 1:Nspecies){ - abGenusSpecies[k,1] ~ dnorm(abTrue[1],tauGenusSpecies[1]) - abGenusSpecies[k,2] ~ dnorm(abTrue[2],tauGenusSpecies[2]) - } - - ### Observation model - ## Errors - # given the data and the priors, establish distributions - # for the observation errors sigmaObs - - sigmaObslog10a ~ dgamma( SD_rObs_log10a, SD_muObs_log10a) - sigmaObsb ~ dgamma( SD_rObs_b, SD_muObs_b) - - # create inverse covariance matrix, with negative parameter correlation roObs - roObs ~ dunif(-0.99,0) - CovObs[1,1] <- pow(sigmaObslog10a,2) - CovObs[2,2] <- pow(sigmaObsb,2) - CovObs[1,2] <- roObs * sigmaObslog10a * sigmaObsb - CovObs[2,1] <- CovObs[1,2] - TauObs[1:2,1:2] <- inverse(CovObs[1:2,1:2]) - - ## likelihood - # given the data, the priors and the covariance, - # create multivariate likelihood distributions for log10(a) and b - - for(i in 1:N){ - TauObsI[i,1:2,1:2] <- TauObs[1:2,1:2] * pow(Weights[i],2) # weighted precision - ab[i,1:2] ~ dmnorm(abGenusSpecies[GenusSpecies[i],1:2],TauObsI[i,1:2,1:2]) - } -} - " - - # Write JAGS model - cat(Model, file=JAGSFILE) - # JAGS settings - Nchains = 3 # number of MCMC chains to be used in JAGS - Nburnin = 1e4 # number of burn-in iterations, to be discarded; 1e4 = 10000 iterations for burn-in - Niter = 3e4 # number of iterations after burn-in; 3e4 = 30000 iterations - Nthin = 1e1 # subset of iterations to be used for analysis; 1e1 = every 10th iteration - - # Run JAGS: define data to be passed on in DataJags; - # determine parameters to be returned in Param2Save; - # call JAGS with function Jags() - DataJags = list(ab=cbind(log10(a),b), N=length(a), Weights=wts, Nspecies=Nspecies, GenusSpecies=GenusSpecies, - prior_mean_b=prior_mean_b, prior_tau_b=prior_tau_b, - prior_mean_log10a=prior_mean_log10a, prior_tau_log10a=prior_tau_log10a, - SD_rObs_log10a=SD_rObs_log10a, SD_muObs_log10a=SD_muObs_log10a, - SD_rObs_b=SD_rObs_b, SD_muObs_b=SD_muObs_b, - SD_rGS_log10a=SD_rGS_log10a, SD_muGS_log10a=SD_muGS_log10a, - SD_rGS_b=SD_rGS_b, SD_muGS_b=SD_muGS_b) - Params2Save = c("abTrue","abGenusSpecies","sigmaGSlog10a","sigmaGSb","sigmaObslog10a","sigmaObsb","roObs") - Jags <- jags(inits=NULL, model.file=JAGSFILE, working.directory=NULL, data=DataJags, - parameters.to.save=Params2Save, n.chains=Nchains, n.thin=Nthin, n.iter=Niter, n.burnin=Nburnin) - Jags$BUGSoutput # contains the results from the JAGS run - - # Analyze output for the relatives - abTrue <- Jags$BUGSoutput$sims.list$abTrue - R_mean_log10a <- mean(abTrue[,1]) # true mean of log10(a) - R_sd_log10a <- sd(abTrue[,1]) # true SE of log10(a) - R_mean_b <- mean(abTrue[,2]) # true mean of b - R_sd_b <- sd(abTrue[,2]) # true SE of b - - # Analyze output for the target species - abGenusSpecies <- Jags$BUGSoutput$sims.list$abGenusSpecies - mean_log10a <- mean(abGenusSpecies[,1,1]) # true mean of log10(a) for the first species= target species - sd_log10a <- sd(abGenusSpecies[,1,1]) # true SE of log10(a) - mean_b <- mean(abGenusSpecies[,1,2]) # true mean of b - sd_b <- sd(abGenusSpecies[,1,2]) # true SE of b - mean_sigma_log10a <- mean(Jags$BUGSoutput$sims.list$sigmaObslog10a) # measurement error of log10(a) - sd_sigma_log10a <- apply(as.matrix(Jags$BUGSoutput$sims.list$sigmaObslog10a), 2, sd) - mean_sigma_b <- mean(Jags$BUGSoutput$sims.list$sigmaObsb) # measurement error of b - sd_sigma_b <- apply(as.matrix(Jags$BUGSoutput$sims.list$sigmaObsb), 2, sd) - ro_ab <- mean(Jags$BUGSoutput$sims.list$roObs) # measurement correlation of log10(a),b - - out.list <- list(N=length(a), mean_log10a=mean_log10a, sd_log10a=sd_log10a, mean_b=mean_b, sd_b=sd_b, - R_mean_log10a=R_mean_log10a, R_sd_log10a=R_sd_log10a, R_mean_b=R_mean_b, R_sd_b=R_sd_b) - return(out.list) - } - -#----------------------------------------------------------------------------- -# Function to do a Bayesian LWR analysis with studies for target species only -#----------------------------------------------------------------------------- -SpecLWR <- function(a, b, wts, prior_mean_b, prior_tau_b, - prior_mean_log10a, prior_tau_log10a, SD_rObs_log10a, SD_muObs_log10a, - SD_rObs_b, SD_muObs_b, SD_rGS_log10a, SD_muGS_log10a, - SD_rGS_b, SD_muGS_b){ - - # Define JAGS model - Model = " - model { - sigma1 ~ dgamma( SD_rObs_log10a, SD_muObs_log10a) # posterior distribution for measurement error in log10a - sigma2 ~ dgamma( SD_rObs_b, SD_muObs_b) # posterior distribution for measurement error in log10a - - ro ~ dunif(-0.99,0) # uniform prior for negative correlation between log10a and b - abTrue[1] ~ dnorm(prior_mean_log10a,prior_tau_log10a) # normal posterior distribution for log10a - abTrue[2] ~ dnorm(prior_mean_b,prior_tau_b) # normal posterior distribution for b - CovObs[1,1] <- pow(sigma1,2) - CovObs[2,2] <- pow(sigma2,2) - CovObs[1,2] <- ro * sigma1 * sigma2 - CovObs[2,1] <- CovObs[1,2] - TauObs[1:2,1:2] <- inverse(CovObs[1:2,1:2]) # create inverse covariance matrix - for(i in 1:N){ - TauObsI[i,1:2,1:2] <- TauObs[1:2,1:2] * pow(Weights[i],2) # converts prior SD into prior weighted precision - - # given the data, the priors and the covariance, create multivariate normal posteriors for log(a) and b - ab[i,1:2] ~ dmnorm(abTrue[1:2],TauObsI[i,1:2,1:2]) - } - } - " - -# Write JAGS model -cat(Model, file=JAGSFILE) -# JAGS settings -Nchains = 3 # number of MCMC chains to be used in JAGS -Nburnin = 1e4 # number of burn-in runs, to be discarded; 10000 iterations for burn-in -Niter = 3e4 # number of iterations after burn-in; 3e4 = 30000 iterations -Nthin = 1e1 # subset of iterations to be used for analysis; 1e1 = every 10th iteration -# Run JAGS: define data to be passed on in DataJags; determine parameters to be returned in Param2Save; call JAGS with function Jags() -DataJags = list(ab=cbind(log10(a),b), N=length(a), Weights=wts, prior_mean_b=prior_mean_b, - prior_tau_b=prior_tau_b, prior_mean_log10a=prior_mean_log10a, prior_tau_log10a=prior_tau_log10a, - SD_rObs_log10a=SD_rObs_log10a, SD_muObs_log10a=SD_muObs_log10a, - SD_rObs_b=SD_rObs_b, SD_muObs_b=SD_muObs_b) -Params2Save = c("abTrue","sigma1","sigma2","ro") -Jags <- jags(inits=NULL, model.file=JAGSFILE, working.directory=NULL, data=DataJags, parameters.to.save=Params2Save, n.chains=Nchains, n.thin=Nthin, n.iter=Niter, n.burnin=Nburnin) -Jags$BUGSoutput # contains the results from the JAGS run -# Analyze output -abTrue <- Jags$BUGSoutput$sims.list$abTrue -mean_log10a <- mean(abTrue[,1]) # true mean of log10(a) -sd_log10a <- sd(abTrue[,1]) # true SE of log10(a) -mean_b <- mean(abTrue[,2]) # true mean of b -sd_b <- sd(abTrue[,2]) # true SE of b -mean_sigma_log10a <- mean(Jags$BUGSoutput$sims.list$sigma1) # measurement error of log10(a) -sd_sigma_log10a <- apply(as.matrix(Jags$BUGSoutput$sims.list$sigma1), 2, sd) -mean_sigma_b <- mean(Jags$BUGSoutput$sims.list$sigma2) # measurement error of b -sd_sigma_b <- apply(as.matrix(Jags$BUGSoutput$sims.list$sigma2), 2, sd) -ro_ab <- mean(Jags$BUGSoutput$sims.list$ro) # measurement correlation of log10(a),b - -out.list <- list(N=length(a), mean_log10a=mean_log10a, sd_log10a=sd_log10a, mean_b=mean_b, sd_b=sd_b) -return(out.list) - -} # End of Functions section - -#-------------------------------- -# Analysis by Family -#-------------------------------- -# Do LWR analysis by Family, Subfamily and Body shape, depending on available LWR studies -# for(Fam in "Acanthuridae") { -for(Fam in Families) { - Subfamilies <- sort(unique(Data$Subfamily[Data$Family==Fam])) - for(SF in Subfamilies) { - for(BS in Bshape) { - # get species (SpecCodes) in this Subfamily and with this body shape - SpecCode.SF.BS <- unique(Data$SpecCode[Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS]) - # if there are species with this body shape - if(length(SpecCode.SF.BS > 0)) { - # get priors for this body shape - prior <- Get.BS.pr(BS) - # get LWR data for this body shape - b_raw <- Data$b[Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS] - cat("\n") - cat("Family =", Fam, ", Subfamily =", SF, ", Body shape =", BS, ", Species =", length(SpecCode.SF.BS), ", LWR =", - length(b_raw[is.na(b_raw)==F]), "\n") - # if no LWR studies exist for this body shape, assign the respective priors to all species - if(length(b_raw[is.na(b_raw)==F])==0) { - # assign priors to species with no LWR in this Subfamily with this body shape - cat("Assigning overall body shape prior to", length(SpecCode.SF.BS), " species \n") - for(SpC in SpecCode.SF.BS) { - out.prior <- data.frame(Fam, SF, BS, SpC, 0, prior$mean_log10a, prior$sd_log10a, prior$mean_b, prior$sd_b, - paste("all LWR estimates for this BS")) - write.table(out.prior, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F) - } - } else { - - # Update priors for this body shape using existing LWR studies - # get LWR data for this Subfamily and body shape - Keep <- which(Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F & Data$Score>0) - wts <- Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses) - a <- Data$a[Keep] - b <- Data$b[Keep] - GenSpec <- paste(Data$Genus[Keep],Data$Species[Keep]) - # add a first dummy record with prior LWR and low score = 0.3, as pseudo target species - # Name of dummy target species is Dum1 dum1 - TargetSpec = paste("Dum1", "dum1") - wts <- c(0.3, wts) - a <- c(10^(prior$mean_log10a), a) - b <- c(prior$mean_b, b) - GenSpec <- c(TargetSpec, GenSpec) - # Relabel GenSpec so that TargetSpec = level 1 - OtherSpecies = unique(GenSpec[GenSpec != TargetSpec]) - GenusSpecies = factor(GenSpec, levels=c(TargetSpec, OtherSpecies)) - Nspecies = nlevels(GenusSpecies) # number of species - # run Bayesian analysis for pseudo target species with Subfamily members - # The resulting R_mean_log10a, R_sd_log10a, R_mean_b, R_sd_b will be used for species without LWR - cat("Updating Subfamily-Bodyshape prior using", Nspecies-1, "species with LWR studies \n") - prior.SFam.BS <- SpecRelLWR(a, b, wts, GenusSpecies, Nspecies, prior_mean_b=prior$mean_b, - prior_tau_b=prior$tau_b, prior_mean_log10a=prior$mean_log10a, - prior_tau_log10a=prior$tau_log10a, SD_rObs_log10a=prior$SD_rObs_log10a, - SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b, - SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a, - SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b, - SD_muGS_b=prior$SD_muGS_b) - - #------------------------------------------------------------------------------------------ - # if there are Genera with >= 5 species with LWR, update body shape priors for these Genera - #------------------------------------------------------------------------------------------ - Genera <- unique(as.character(Data$Genus[Keep])) - # create empty list of lists for storage of generic priors - prior.Gen.BS <- rep(list(list()),length(Genera)) # create a list of empty lists - names(prior.Gen.BS) <- Genera # name the list elements according to the Genera - for(Genus in Genera){ - # check if Genus contains >= 5 species with LWR data - if(length(unique(Data$SpecCode[Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F & - Data$Score>0 & Data$Genus==Genus]))>=5) { - # run Subfamily analysis with only data for this genus - Keep <- which(Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F & Data$Score>0 & - Data$Genus==Genus) - wts <- Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses) - a <- Data$a[Keep] - b <- Data$b[Keep] - GenSpec <- paste(Data$Genus[Keep],Data$Species[Keep]) - # add a first dummy record with prior LWR and low score = 0.3, as pseudo target species - # Name of dummy target species is Dum1 dum1 - TargetSpec = paste("Dum1", "dum1") - wts <- c(0.3, wts) - a <- c(10^(prior$mean_log10a), a) - b <- c(prior$mean_b, b) - GenSpec <- c(TargetSpec, GenSpec) - # Relabel GenSpec so that TargetSpec = level 1 - OtherSpecies = unique(GenSpec[GenSpec != TargetSpec]) - GenusSpecies = factor(GenSpec, levels=c(TargetSpec, OtherSpecies)) - Nspecies = nlevels(GenusSpecies) # number of species - # run Bayesian analysis for pseudo target species with Genus members - # R_mean_log10a, R_sd_log10a, R_mean_b, R_sd_b will be used for species without LWR - cat("Updating prior for Genus =", Genus, ", with", Nspecies -1, "LWR Species \n") - prior.Gen.BS[[Genus]] <- SpecRelLWR(a, b, wts, GenusSpecies, Nspecies, - prior_mean_b=prior.SFam.BS$R_mean_b, - prior_tau_b=1/prior.SFam.BS$R_sd_b^2, - prior_mean_log10a=prior.SFam.BS$R_mean_log10a, - prior_tau_log10a=1/prior.SFam.BS$R_sd_log10a, - SD_rObs_log10a=prior$SD_rObs_log10a, - SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b, - SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a, - SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b, - SD_muGS_b=prior$SD_muGS_b) - } - } - # new Subfamily-BS priors have been generated - # for some genera, new Genus-BS priors have been generated - # --------------------------------------------------------------------- - # Loop through all species in this Subfamily-BS; assign LWR as appropriate - # --------------------------------------------------------------------- - for(SpC in SpecCode.SF.BS) { - Genus <- as.character(unique(Data$Genus[Data$SpecCode==SpC])) - Species <- as.character(unique(Data$Species[Data$SpecCode==SpC])) - TargetSpec = paste(Genus, Species) - LWR <- length(Data$b[Data$SpecCode==SpC & is.na(Data$b)==F & Data$Score>0]) - LWRGenspec <- length(unique(Data$SpecCode[Data$BodyShapeI==BS & is.na(Data$b)==F & - Data$Score>0 & Data$Genus==Genus])) - LWRSFamspec <- length(unique(Data$SpecCode[Data$BodyShapeI==BS & is.na(Data$b)==F & - Data$Score>0 & Data$Family==Fam & Data$Subfamily==SF])) - #--------------------------------------------------------- - # >= 5 LWR in target species, run single species analysis - if(LWR >= 5) { - # Run analysis with data only for this species - Keep <- which(Data$SpecCode==SpC & is.na(Data$b)==F & Data$Score>0) - wts = Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses) - a = Data$a[Keep] - b = Data$b[Keep] - - # determine priors to be used - if(LWRGenspec >= 5) { - prior_mean_b=prior.Gen.BS[[Genus]]$R_mean_b - prior_tau_b=1/prior.Gen.BS[[Genus]]$R_sd_b^2 - prior_mean_log10a=prior.Gen.BS[[Genus]]$R_mean_log10a - prior_tau_log10a=1/prior.Gen.BS[[Genus]]$R_sd_log10a^2 - } else - if (LWRSFamspec > 0) { - prior_mean_b=prior.SFam.BS$R_mean_b - prior_tau_b=1/prior.SFam.BS$R_sd_b^2 - prior_mean_log10a=prior.SFam.BS$R_mean_log10a - prior_tau_log10a=1/prior.SFam.BS$R_sd_log10a^2 - } else { - prior_mean_b=prior$mean_b - prior_tau_b=prior$tau_b - prior_mean_log10a=prior$mean_log10a - prior_tau_log10a=prior$tau_log10a - } - cat("Running single species analysis for", TargetSpec, "LWR =", LWR, ", LWR species in Genus=",LWRGenspec,"\n" ) - # call function for single species analysis - post <- SpecLWR(a, b, wts, prior_mean_b=prior_mean_b, - prior_tau_b=prior_tau_b, prior_mean_log10a=prior_mean_log10a, - prior_tau_log10a=prior_tau_log10a, SD_rObs_log10a=prior$SD_rObs_log10a, - SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b, - SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a, - SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b, - SD_muGS_b=prior$SD_muGS_b) - out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(post$mean_log10a, digits=3), format(post$sd_log10a, digits=3), format(post$mean_b, disgits=3), format(post$sd_b, digits=3), - paste("LWR estimates for this species")) - write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F) - - } else - #-------------------------------------------------------- - # 1-4 LWR in target species and >= 5 LWR species in Genus - # run hierarchical analysis for genus members, with Subfamily-BS prior - if(LWR >= 1 & LWRGenspec >=5) { - # run Subfamily analysis with only data for this genus - Keep <- which(Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F & Data$Score>0 & - Data$Genus==Genus) - wts <- Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses) - a <- Data$a[Keep] - b <- Data$b[Keep] - GenSpec <- paste(Data$Genus[Keep],Data$Species[Keep]) - - # Relabel GenSpec so that TargetSpec = level 1 - OtherSpecies = unique(GenSpec[GenSpec != TargetSpec]) - GenusSpecies = factor(GenSpec, levels=c(TargetSpec, OtherSpecies)) - Nspecies = nlevels(GenusSpecies) # number of species - # run Bayesian analysis for target species with Genus members - cat("Running analysis with congeners for", TargetSpec, ", LWR =", LWR,", LWR species in Genus =", LWRGenspec,"\n") - post <- SpecRelLWR(a, b, wts, GenusSpecies, Nspecies, - prior_mean_b=prior.SFam.BS$R_mean_b, - prior_tau_b=1/prior.SFam.BS$R_sd_b^2, - prior_mean_log10a=prior.SFam.BS$R_mean_log10a, - prior_tau_log10a=1/prior.SFam.BS$R_sd_log10a^2, - SD_rObs_log10a=prior$SD_rObs_log10a, - SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b, - SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a, - SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b, - SD_muGS_b=prior$SD_muGS_b) - out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(post$mean_log10a, digits=3), format(post$sd_log10a, digits=3), format(post$mean_b, disgits=3), format(post$sd_b, digits=3), - paste("LWR estimates for species & Genus-BS")) - write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F) - } else - - #------------------------------------------------------- - # 1-4 LWR in target species and < 5 LWR species in Genus - # run hierarchical analysis for Subfamily members, with bodyshape prior - - if(LWR >= 1 & LWRSFamspec > 1) { - # run Subfamily analysis - Keep <- which(Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F & Data$Score>0) - wts <- Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses) - a <- Data$a[Keep] - b <- Data$b[Keep] - GenSpec <- paste(Data$Genus[Keep],Data$Species[Keep]) - # Relabel GenSpec so that TargetSpec = level 1 - OtherSpecies = unique(GenSpec[GenSpec != TargetSpec]) - GenusSpecies = factor(GenSpec, levels=c(TargetSpec, OtherSpecies)) - Nspecies = nlevels(GenusSpecies) # number of species - # run Bayesian analysis for target species with Subfamily members - cat("Running analysis with Subfamily members for", TargetSpec, ", LWR =", LWR,", LWR species in Subfamily-BS =", - LWRSFamspec, "\n") - post <- SpecRelLWR(a, b, wts, GenusSpecies, Nspecies, - prior_mean_b=prior$mean_b, - prior_tau_b=prior$tau_b, - prior_mean_log10a=prior$mean_log10a, - prior_tau_log10a=prior$tau_log10a, - SD_rObs_log10a=prior$SD_rObs_log10a, - SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b, - SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a, - SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b, - SD_muGS_b=prior$SD_muGS_b) - out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(post$mean_log10a, digits=3), format(post$sd_log10a, digits=3), - format(post$mean_b, disgits=3), format(post$sd_b, digits=3), - paste("LWR estimates for species & Subfamily-BS")) - write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F) - } else - #-------------------------------------------------- - # assign Genus-BS priors to target species - if(LWRGenspec >= 5) { - cat("Assign Genus-BS prior for", TargetSpec, "\n") - out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(prior.Gen.BS[[Genus]]$mean_log10a, digits=3), - format(prior.Gen.BS[[Genus]]$sd_log10a, digits=3), - format(prior.Gen.BS[[Genus]]$mean_b, digits=3), format(prior.Gen.BS[[Genus]]$sd_b, digits=3), - paste("LWR estimates for this Genus-BS")) - write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F) - } else { - # ----------------------------------------------- - # assign Subfamily-BS priors to target species - cat("Assign Subfamily-BS prior for", TargetSpec,"\n") - out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(prior.SFam.BS$mean_log10a, digits=3), format(prior.SFam.BS$sd_log10a, digits=3), - format(prior.SFam.BS$mean_b, digits=3), format(prior.SFam.BS$sd_b, digits=3), paste("LWR estimates for this Subfamily-BS")) - write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F) - } - } # end of species loop for this Subfamily and body shape - - } # end of section dealing with Subfamily - body shapes that contain LWR estimates - - } # end of section that deals with Subfamily - body shapes that contain species - - } # end of body shape section - - } # end of Subfamily section - -} # end of Family section -cat("End", date(),"\n") - - - - - - - - - diff --git a/PARALLEL_PROCESSING/YASMEEN-matcher-1.2.0.1.jar b/PARALLEL_PROCESSING/YASMEEN-matcher-1.2.0.1.jar deleted file mode 100644 index ef14b6f..0000000 Binary files a/PARALLEL_PROCESSING/YASMEEN-matcher-1.2.0.1.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/YASMEEN-parser-1.2.0.jar b/PARALLEL_PROCESSING/YASMEEN-parser-1.2.0.jar deleted file mode 100644 index 42e28c2..0000000 Binary files a/PARALLEL_PROCESSING/YASMEEN-parser-1.2.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/c3p0-0.9.1.2.jar b/PARALLEL_PROCESSING/c3p0-0.9.1.2.jar deleted file mode 100644 index 0f42d60..0000000 Binary files a/PARALLEL_PROCESSING/c3p0-0.9.1.2.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/common-configuration-scanner-1.0.1-SNAPSHOT.jar b/PARALLEL_PROCESSING/common-configuration-scanner-1.0.1-SNAPSHOT.jar deleted file mode 100644 index d5c087d..0000000 Binary files a/PARALLEL_PROCESSING/common-configuration-scanner-1.0.1-SNAPSHOT.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/common-encryption-1.0.1-3.5.0.jar b/PARALLEL_PROCESSING/common-encryption-1.0.1-3.5.0.jar deleted file mode 100644 index 20f4a73..0000000 Binary files a/PARALLEL_PROCESSING/common-encryption-1.0.1-3.5.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/common-gcore-resources-1.2.0-3.5.0.jar b/PARALLEL_PROCESSING/common-gcore-resources-1.2.0-3.5.0.jar deleted file mode 100644 index 27a0ec5..0000000 Binary files a/PARALLEL_PROCESSING/common-gcore-resources-1.2.0-3.5.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/common-gcore-stubs-1.2.0-3.5.0.jar b/PARALLEL_PROCESSING/common-gcore-stubs-1.2.0-3.5.0.jar deleted file mode 100644 index caf569a..0000000 Binary files a/PARALLEL_PROCESSING/common-gcore-stubs-1.2.0-3.5.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/common-scope-1.2.1-SNAPSHOT.jar b/PARALLEL_PROCESSING/common-scope-1.2.1-SNAPSHOT.jar deleted file mode 100644 index 2602d3e..0000000 Binary files a/PARALLEL_PROCESSING/common-scope-1.2.1-SNAPSHOT.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/common-scope-maps-1.0.2-3.5.0.jar b/PARALLEL_PROCESSING/common-scope-maps-1.0.2-3.5.0.jar deleted file mode 100644 index e1ee05d..0000000 Binary files a/PARALLEL_PROCESSING/common-scope-maps-1.0.2-3.5.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/commons-collections-3.1.jar b/PARALLEL_PROCESSING/commons-collections-3.1.jar deleted file mode 100644 index 41e230f..0000000 Binary files a/PARALLEL_PROCESSING/commons-collections-3.1.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/commons-io-1.2.jar b/PARALLEL_PROCESSING/commons-io-1.2.jar deleted file mode 100644 index b2867cd..0000000 Binary files a/PARALLEL_PROCESSING/commons-io-1.2.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/discovery-client-1.0.1-3.5.0.jar b/PARALLEL_PROCESSING/discovery-client-1.0.1-3.5.0.jar deleted file mode 100644 index ea61c43..0000000 Binary files a/PARALLEL_PROCESSING/discovery-client-1.0.1-3.5.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/dom4j-1.6.1.jar b/PARALLEL_PROCESSING/dom4j-1.6.1.jar deleted file mode 100644 index c8c4dbb..0000000 Binary files a/PARALLEL_PROCESSING/dom4j-1.6.1.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/ecological-engine-1.8.6-SNAPSHOT.jar b/PARALLEL_PROCESSING/ecological-engine-1.8.6-SNAPSHOT.jar deleted file mode 100644 index a660bd0..0000000 Binary files a/PARALLEL_PROCESSING/ecological-engine-1.8.6-SNAPSHOT.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/ecological-engine-smart-executor-1.1.0-SNAPSHOT.jar b/PARALLEL_PROCESSING/ecological-engine-smart-executor-1.1.0-SNAPSHOT.jar deleted file mode 100644 index 3d6ddc5..0000000 Binary files a/PARALLEL_PROCESSING/ecological-engine-smart-executor-1.1.0-SNAPSHOT.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/hibernate3.jar b/PARALLEL_PROCESSING/hibernate3.jar deleted file mode 100644 index f97dcdb..0000000 Binary files a/PARALLEL_PROCESSING/hibernate3.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/ic-client-1.0.1-3.5.0.jar b/PARALLEL_PROCESSING/ic-client-1.0.1-3.5.0.jar deleted file mode 100644 index 792b11f..0000000 Binary files a/PARALLEL_PROCESSING/ic-client-1.0.1-3.5.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/jaxen-1.1.2.jar b/PARALLEL_PROCESSING/jaxen-1.1.2.jar deleted file mode 100644 index 69de309..0000000 Binary files a/PARALLEL_PROCESSING/jaxen-1.1.2.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/jta-1.1.jar b/PARALLEL_PROCESSING/jta-1.1.jar deleted file mode 100644 index 6d225b7..0000000 Binary files a/PARALLEL_PROCESSING/jta-1.1.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/log4j-1.2.16.jar b/PARALLEL_PROCESSING/log4j-1.2.16.jar deleted file mode 100644 index 3f9d847..0000000 Binary files a/PARALLEL_PROCESSING/log4j-1.2.16.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/mongo-java-driver-2.12.4.jar b/PARALLEL_PROCESSING/mongo-java-driver-2.12.4.jar deleted file mode 100644 index ecbaecd..0000000 Binary files a/PARALLEL_PROCESSING/mongo-java-driver-2.12.4.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/postgresql-8.4-702.jdbc4.jar b/PARALLEL_PROCESSING/postgresql-8.4-702.jdbc4.jar deleted file mode 100644 index 8b0c761..0000000 Binary files a/PARALLEL_PROCESSING/postgresql-8.4-702.jdbc4.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/slf4j-api-1.6.0.jar b/PARALLEL_PROCESSING/slf4j-api-1.6.0.jar deleted file mode 100644 index db92f9a..0000000 Binary files a/PARALLEL_PROCESSING/slf4j-api-1.6.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/slf4j-log4j12-1.6.0.jar b/PARALLEL_PROCESSING/slf4j-log4j12-1.6.0.jar deleted file mode 100644 index 42cca51..0000000 Binary files a/PARALLEL_PROCESSING/slf4j-log4j12-1.6.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/storage-manager-core-2.2.0-3.7.0.jar b/PARALLEL_PROCESSING/storage-manager-core-2.2.0-3.7.0.jar deleted file mode 100644 index d38e338..0000000 Binary files a/PARALLEL_PROCESSING/storage-manager-core-2.2.0-3.7.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/storage-manager-wrapper-2.2.0-3.7.0.jar b/PARALLEL_PROCESSING/storage-manager-wrapper-2.2.0-3.7.0.jar deleted file mode 100644 index bc38dfb..0000000 Binary files a/PARALLEL_PROCESSING/storage-manager-wrapper-2.2.0-3.7.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/xalan-2.6.0.jar b/PARALLEL_PROCESSING/xalan-2.6.0.jar deleted file mode 100644 index 73cf175..0000000 Binary files a/PARALLEL_PROCESSING/xalan-2.6.0.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/xpp3_min-1.1.4c.jar b/PARALLEL_PROCESSING/xpp3_min-1.1.4c.jar deleted file mode 100644 index 813a9a8..0000000 Binary files a/PARALLEL_PROCESSING/xpp3_min-1.1.4c.jar and /dev/null differ diff --git a/PARALLEL_PROCESSING/xstream-1.3.1.jar b/PARALLEL_PROCESSING/xstream-1.3.1.jar deleted file mode 100644 index 4ef4219..0000000 Binary files a/PARALLEL_PROCESSING/xstream-1.3.1.jar and /dev/null differ