ELEFAN_Boot_c7
Ralf Schwamborn
26 de abril de 2018
R Markdown
# ELEFAN_BOOT (Version 1.0.XXII.c7, beta4)
# Analysis of LFD data: from RSA to confidence ellipses
# Version XXII_beta4
# for LINUX_server 1 - (no plots, no beeps)
# Copyright: Ralf Schwamborn, November 2017
# Only necessary code for LINUX SERVER BOOTSTRAPS (no plots, no beeps)
# Experiments 1,2,3,4
# Experiment 1 (OK): Boot1 vs FULL Boot2, with _GA, alba
# Experiment 2 (OK): Boot1 vs FULL Boot2, with _SA, alba
# Experiment 3 (OK): Boot1 vs FULL Boot2, with _GA, synLF5b (12 months only)
# Experiment 4 (OK): Boot1 vs FULL Boot2, with _SA, synLF5b (12 months only)
# Experiment 5 (not implem.): Boot1 vs FULL Boot2, using a RSA plot, alba
# Experiment 6 (not implem.): Boot1 vs FULL Boot2, using a RSA plot, synLF5b (12 months only) ]
# [Experiment x (not yet ready): effect of seed distrib on output distrib, rnorm vs runif ]
### STEP 1: Download and install TropFishR and other necessary packages
### (this step takes a lot of time in LINUX ..., but is fast in Windows or Apple)
install.packages("parallel", repos='http://cran.us.r-project.org') # OK
install.packages("doParallel", repos='http://cran.us.r-project.org') # OK
install.packages("doRNG", repos='http://cran.us.r-project.org') # OK
install.packages("TropFishR", repos='http://cran.us.r-project.org') # OK, but old version
install.packages("openssl-devel", repos='http://cran.us.r-project.org') # oops, not working, server blocks OpenSSL
install.packages("libssl-dev", repos='http://cran.us.r-project.org') # oops, not working, server blocks OpenSSL
install.packages("openssl ", repos='http://cran.us.r-project.org') # oops, not working, server blocks OpenSSL
install.packages("git2r", repos='http://cran.us.r-project.org') # not working, server blocks OpenSSL
install.packages("devtools", repos='http://cran.us.r-project.org') # oops, not working, server blocks OpenSSL
devtools::install_github("tokami/TropFishR", ref="smallImpros", force=TRUE) # oops, not working, server blocks OpenSSL
### STEP 2: load TropFishR and other necessary packages, prepare a "alba2" file for use, with one additional (empty) large size class
library(TropFishR)
library(parallel)
library(doParallel)
library(doRNG)
# setwd("~/WG_New_Length-Based_Methods")
alba2 <- alba
cat.t <- alba2$catch
colnames(cat.t)
dim(cat.t) # matrix dimensions 14 rows x 7 cols
zeros.mat <- matrix(0, nrow = 2, ncol = 7)
cat.t2 <- rbind(cat.t,zeros.mat )
dim(cat.t2) # new matrix dimensions 16 rows x 7 cols
alba2$sample.no <- c( 1, 2, 3 ,4 ,5 ,6 ,7 ,8, 9, 10, 11, 12, 13 ,14 ,
15, 16)
alba2$midLengths = c(1.5 ,2.5 ,3.5, 4.5, 5.5, 6.5 ,7.5 ,8.5 ,9.5, 10.5, 11.5, 12.5, 13.5 ,14.5,
15.5, 16.5)
alba2$catch <- cat.t2
saveRDS(alba2,file = "alba2.rds") # saves "alba2.rds" to disk
### STEP 3: RUN BOOTSTRAP EXPERIMENTS
#######################################################
# Experiment 1
#
# Boot1 vs FULL Boot2, with _GA, alba
##################
# BOOTSTRAP 1
# 10 runs take about 1 minute
# RUN BOOSTRAP 1 # maxit = 8, 10 runs = 1 MINUTE, sample_seeds <- round(runif(100, max = 100, min = 1),0)
for (i in 1:400) {
dat.time<- Sys.time() # Date and Time
dat.time.num<- as.numeric(dat.time) # Date and Time (numeric)
seed_glob <- i + dat.time.num # Global Seed value (date+time+"i"), used to Reset the Random Number Generator
set.seed(seed_glob) # Resets the Random Number Generator
library(TropFishR)
library(parallel)
library(doParallel)
library(doRNG)
Linf <- rep(NA, 1)
K <-rep(NA, 1)
t_anchor <-rep(NA, 1)
Phi_prime <- rep(NA, 1)
Rn_max <- rep(NA, 1)
seed_fit<- rep(NA, 1)
outp.table <- data.frame(Linf, K, t_anchor, Phi_prime, Rn_max,
dat.time,dat.time.num,i, seed_glob, seed_fit)
alba2 <- readRDS("alba2.rds")
sample_seeds <- round(runif(100, max = 100, min = 1),0)
seed_fit <- sample(sample_seeds,1)
output <- ELEFAN_GA(alba2, seasonalised = F,
low_par = list(Linf = 5, K = 0.1, t_anchor = 0, C = 0, ts= 0),
up_par = list(Linf = 15, K = 7, t_anchor = 1, C = 0, ts = 1),
popSize = 40, maxiter = 8,
MA = 9, plot = F, plot.score = F , seed = seed_fit,
parallel = T, beep = F)
output$par
output$Rn_max
outp.table[1,1] <- round(as.numeric(output$par[1]), 2)
outp.table[1,2] <- round(as.numeric(output$par[2]), 3)
outp.table[1,3] <- round(as.numeric(output$par[3]), 3)
outp.table[1,4] <- round(as.numeric(output$par[4]), 3)
outp.table[1,5] <- round(output$Rn_max, 4)
outp.table[1,10] <- seed_fit
#setwd("~/WG_New_Length-Based_Methods")
write.table(outp.table,file = "exp1BOOT_1_iter08_outp_table_GAcalbal_setSeed_bigtable.csv" ,
append = TRUE,sep= ";", col.names = FALSE, row.names = FALSE)
# clear memory (avoids crashing the R session)
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
}
# setwd("~/WG_New_Length-Based_Methods")
boot.out11GA <- read.table(file = "exp1BOOT_1_iter08_outp_table_GAcalbal_setSeed_bigtable.csv" , sep= ";")
names(boot.out11GA) <-c("Linf",
"K",
"t_anchor",
"Phi_p" ,
"Rn_max",
"dat.time", "dat.time.num", "i", "seed_glob",
"seed_fit")
# Best fit parameters
boot.out11GA[which.max(boot.out11GA$Rn_max),]
# Conf.ints
quantile(boot.out11GA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out11GA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out11GA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
length(boot.out11GA$Linf) # n = 919 successful runs
plot(boot.out11GA$Linf)
plot(boot.out11GA$K)
plot(boot.out11GA$Phi_p)
##################
# FULL BOOTSTRAP 2
#
# 10 runs take about 1 minute
for (i in 1:400)
{
dat.time<- Sys.time() # Date and Time
dat.time.num<- as.numeric(dat.time) # Date and Time (numeric)
seed_glob <- i + dat.time.num # Global Seed value (date+time+"i"), used to Reset the Random Number Generator
set.seed(seed_glob) # Resets the Random Number Generator
# Sample and replace ALL monthly LFQ data
library(TropFishR)
library(parallel)
library(doParallel)
library(doRNG)
K.seeds_alb = round(rnorm(100, 1.4, sd = 0.05),2)
Linf.seeds_alb = round(rnorm(100, 10, sd = 0.5),1)
alba2 <- readRDS("alba2.rds")
alba2.t <- alba2
ncol(alba2.t$catch) # n = x montlhy LFQ data
mlen <- alba2.t$midLengths # midlengths
c.int <- (mlen[2]-mlen[1]) #class interval (1 cm)
# "m" ist the month(i.e, the catch file column), 1 to 25
for (m in 1:(ncol(alba2.t$catch))) {
catch.m <- alba2.t$catch[,m]
#plot(catch.m ~ mlen, xlim = c(20,90),main = "original LFD")# original LFD
#rebuild original data
reb.catch.m <- rep(mlen,catch.m) # recreates a proxy of original data
# take a sample
sam.reb.catch.m<- sample(reb.catch.m,replace = TRUE)
hist.sam.reb.catch.m <- hist(sam.reb.catch.m, breaks = c(mlen-(c.int/2)), plot = FALSE)
#replaces data in LFQ file "2" with bootstrap sample data taken from "2"
alba2.t$catch[1:length(hist.sam.reb.catch.m$counts) ,m] <- hist.sam.reb.catch.m$counts
}
alba2.run1 <- alba2.t # new sample
sample_seeds <- round(runif(100, max = 100, min = 1),0)
output2.new <- ELEFAN_GA(alba2.run1, seasonalised = F,
low_par = list(Linf = 5, K = 0.1, t_anchor = 0, C = 0, ts= 0),
up_par = list(Linf = 15, K = 7, t_anchor = 1, C = 0, ts = 1),
popSize = 40, maxiter = 8,
MA = 9, plot = F, plot.score = F ,seed = sample(sample_seeds,1),
parallel = T , beep = F)
output2.new$par
output2.new$Rn_max
output2.new$par
#plot(output2.new)
#date.time <- date()
Linf_out <- round(as.numeric(output2.new$par[1]), 2)
K_out <- round(as.numeric(output2.new$par[2]), 3)
t_anchor <- round(as.numeric(output2.new$par[3]), 3)
Phi_p <- round(as.numeric(output2.new$par[4]), 3)
Rn_max <- round(as.numeric(output2.new$Rn_max),3)
output.table <- data.frame(Linf_out, K_out, t_anchor, Phi_p,Rn_max)
write.table(output.table, file = "exp1BOOT_2_iter08_outp_table_GAcalbah.csv",
sep = ";",col.names = FALSE, row.names = FALSE, append = TRUE )
rm(list=ls())
gc()
} # end of loop
boot.out12GA <- read.table( file = "exp1BOOT_2_iter08_outp_table_GAcalbah.csv", sep = ";")
names(boot.out12GA) <- c("Linf", "K", "t_anchor", "Phi_p","Rn_max")
# Best fit parameters
boot.out12GA[which.max(boot.out12GA$Rn_max),]
# Conf.ints
quantile(boot.out12GA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out12GA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out12GA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
#######################################################
# Experiment 2
#
# Boot1 vs FULL Boot2, with _SA , alba
##################
# BOOTSTRAP 1
# 10 runs take about 1 minute
# RUN BOOSTRAP 1 # maxit = 8, 10 runs = 1 MINUTE, sample_seeds <- round(runif(100, max = 100, min = 1),0)
for (i in 1:800) {
dat.time<- Sys.time() # Date and Time
dat.time.num<- as.numeric(dat.time) # Date and Time (numeric)
seed_glob <- i + dat.time.num # Global Seed value (date+time+"i"), used to Reset the Random Number Generator
set.seed(seed_glob) # Resets the Random Number Generator
library(TropFishR)
K.seeds_alb = round(runif(100, min = 0.3, max = 4),2)
# ASSIGN a series of Linf seed values for the ELEFAN_GA runs
Linf.seeds_alb = round(runif(100, min = 8,max = 12 ),1)
Linf.out.vec <- rep(NA, 1)
K.out.vec <-rep(NA, 1)
t_anchor.out.vec <-rep(NA, 1)
phi.out.vec <- rep(NA, 1)
Rn_max.score.vec <- rep(NA, 1)
outp.table <- data.frame(Linf.out.vec, K.out.vec, t_anchor.out.vec, phi.out.vec, Rn_max.score.vec)
alba2 <- readRDS("alba2.rds")
sample_seeds <- round(runif(100, max = 100, min = 1),0)
output <- ELEFAN_SA(alba2, seasonalised = F,
init_par = list(Linf = sample(Linf.seeds_alb,1), K = sample(K.seeds_alb,1), t_anchor = 0.5, C = 0, ts = 0),
low_par = list(Linf = 5, K = 0.1, t_anchor = 0, C = 0, ts= 0),
up_par = list(Linf = 15, K = 7, t_anchor = 1, C = 0, ts = 1),
beep = F, maxit = 8,
MA = 9, plot = F, plot.score = F
) # SA_time = 0.5
output$par
output$Rn_max
outp.table[1,1] <- round(as.numeric(output$par[1]), 2)
outp.table[1,2] <- round(as.numeric(output$par[2]), 3)
outp.table[1,3] <- round(as.numeric(output$par[3]), 3)
outp.table[1,4] <- round(as.numeric(output$par[4]), 3)
outp.table[1,5] <- round(output$Rn_max, 4)
#setwd("~/WG_New_Length-Based_Methods")
write.table(outp.table,file = "exp2BOOT_1_iter08_outp_table_SAcalbah.csv" ,
append = TRUE,sep= ";", col.names = FALSE, row.names = FALSE)
# clear memory (avoids crashing the R session)
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
}
# setwd("~/WG_New_Length-Based_Methods")
boot.out21SA <- read.table(file = "exp2BOOT_1_iter08_outp_table_SAcalbah.csv" , sep= ";")
names(boot.out21SA) <-c("Linf",
"K",
"t_anchor",
"Phi_p" ,
"Rn_max" )
# Best fit parameters
boot.out21SA[which.max(boot.out21SA$Rn_max),]
# Conf.ints
quantile(boot.out21SA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out21SA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out21SA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
##################
# FULL BOOTSTRAP 2
#
# 10 runs take about 1 minute
for (i in 1:200)
{
dat.time<- Sys.time() # Date and Time
dat.time.num<- as.numeric(dat.time) # Date and Time (numeric)
seed_glob <- i + dat.time.num # Global Seed value (date+time+"i"), used to Reset the Random Number Generator
set.seed(seed_glob) # Resets the Random Number Generator
# Sample and replace ALL monthly LFQ data
library(TropFishR)
K.seeds_alb = round(rnorm(100, 1.4, sd = 0.05),2)
Linf.seeds_alb = round(rnorm(100, 10, sd = 0.5),1)
alba2 <- readRDS("alba2.rds")
alba2.t <- alba2
ncol(alba2.t$catch) # n = x montlhy LFQ data
mlen <- alba2.t$midLengths # midlengths
c.int <- (mlen[2]-mlen[1]) #class interval (1 cm)
# "m" ist the month(i.e, the catch file column), 1 to 25
for (m in 1:(ncol(alba2.t$catch))) {
catch.m <- alba2.t$catch[,m]
#plot(catch.m ~ mlen, xlim = c(20,90),main = "original LFD")# original LFD
#rebuild original data
reb.catch.m <- rep(mlen,catch.m) # recreates a proxy of original data
# take a sample
sam.reb.catch.m<- sample(reb.catch.m,replace = TRUE)
hist.sam.reb.catch.m <- hist(sam.reb.catch.m, breaks = c(mlen-(c.int/2)), plot = FALSE)
#replaces data in LFQ file "2" with bootstrap sample data taken from "2"
alba2.t$catch[1:length(hist.sam.reb.catch.m$counts) ,m] <- hist.sam.reb.catch.m$counts
}
alba2.run1 <- alba2.t # new sample
output2.new <- ELEFAN_SA(alba2.run1, seasonalised = F,
init_par = list(Linf = sample(Linf.seeds_alb,1), K = sample(K.seeds_alb,1), t_anchor = 0.5, C = 0, ts = 0),
low_par = list(Linf = 5, K = 0.1, t_anchor = 0, C = 0, ts= 0),
up_par = list(Linf = 15, K = 7, t_anchor = 1, C = 0, ts = 1),
beep = F, maxit = 8,
MA = 9, plot = F, plot.score = F
) # SA_time = 0.5
output2.new$par
output2.new$Rn_max
output2.new$par
#plot(output2.new)
#date.time <- date()
Linf_out <- round(as.numeric(output2.new$par[1]), 2)
K_out <- round(as.numeric(output2.new$par[2]), 3)
t_anchor <- round(as.numeric(output2.new$par[3]), 3)
Phi_p <- round(as.numeric(output2.new$par[4]), 3)
Rn_max <- round(as.numeric(output2.new$Rn_max),3)
output.table <- data.frame(Linf_out, K_out, t_anchor, Phi_p,Rn_max )
write.table(output.table, file = "exp2BOOT_2_iter08_outp_table_SAdalbah.csv",
sep = ";",col.names = FALSE, row.names = F,
append = TRUE )
rm(list=ls())
gc()
} # end of loop
# setwd("~/WG_New_Length-Based_Methods")
boot.out22SA <- read.table( file = "exp2BOOT_2_iter08_outp_table_SAdalbah.csv", sep = ";")
names(boot.out22SA) <- c("Linf", "K", "t_anchor", "Phi_p","Rn_max")
# Best fit parameters
boot.out22SA[which.max(boot.out22SA$Rn_max),]
# Conf.ints
quantile(boot.out22SA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out22SA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out22SA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
#######################################################
# Experiment 3
#
# Boot1 vs FULL Boot2, with _GA, Synlf5b
library(TropFishR)
# synLFQ5 # Linf = 80, K= 0.5,t_anchor = 0.25, C = 0, 2 years, perfect data
# use one year only
synLF5b <- synLFQ5
synLF5b$catch <- synLF5b$catch[,1:12]
synLF5b$dates <- synLF5b$dates[1:12]
##################
# BOOTSTRAP 1
# 5 runs take about 1 minute. 10 to 15 secs (averge: 12 secs) per run with two intel cores
# RUN BOOSTRAP 1 # maxit = 20, 10 runs = 1 MINUTE, sample_seeds <- round(runif(100, max = 100, min = 1),0)
for (i in 1:200) {
dat.time<- Sys.time() # Date and Time
dat.time.num<- as.numeric(dat.time) # Date and Time (numeric)
seed_glob <- i + dat.time.num # Global Seed value (date+time+"i"), used to Reset the Random Number Generator
set.seed(seed_glob) # Resets the Random Number Generator
library(TropFishR)
library(parallel)
library(doParallel)
library(doRNG)
synLF5b <- synLFQ5
synLF5b$catch <- synLF5b$catch[,1:12]
synLF5b$dates <- synLF5b$dates[1:12]
Linf <- rep(NA, 1)
K <-rep(NA, 1)
t_anchor <-rep(NA, 1)
Phi_prime <- rep(NA, 1)
Rn_max <- rep(NA, 1)
seed_fit<- rep(NA, 1)
outp.table <- data.frame(Linf, K, t_anchor, Phi_prime, Rn_max,
dat.time,dat.time.num,i, seed_glob, seed_fit)
sample_seeds <- round(rnorm(1000, mean = 200, sd = 10),0)
output <- ELEFAN_GA(synLF5b, seasonalised = F,
low_par = list(Linf = 40, K = 0.1, t_anchor = 0, C = 0, ts= 0),
up_par = list(Linf = 160, K = 7, t_anchor = 1, C = 0, ts = 1),
popSize = 50, maxiter = 20,
MA = 11, plot = F, plot.score = F , seed = sample(sample_seeds,1),
parallel = T, beep = F)
output$par
output$Rn_max
outp.table[1,1] <- round(as.numeric(output$par[1]), 2)
outp.table[1,2] <- round(as.numeric(output$par[2]), 3)
outp.table[1,3] <- round(as.numeric(output$par[3]), 3)
outp.table[1,4] <- round(as.numeric(output$par[4]), 3)
outp.table[1,5] <- round(output$Rn_max, 4)
outp.table[1,10] <- seed_fit
# setwd("~/WG_New_Length-Based_Methods")
write.table(outp.table,file = "exp3BOOT_1_iter20_outp_table_GAd_syn5bmsetseed_bigtable.csv" ,
append = TRUE,sep= ";", col.names = FALSE, row.names = FALSE)
# clear memory (avoids crashing the R session)
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
}
# setwd("~/WG_New_Length-Based_Methods")
boot.out31GA <- read.table(file = "exp3BOOT_1_iter20_outp_table_GAd_syn5bmsetseed_bigtable.csv" , sep= ";")
names(boot.out31GA) <-c("Linf", "K",
"t_anchor",
"Phi_p" ,
"Rn_max" ,
"dat.time", "dat.time.num", "i", "seed_glob",
"seed_fit")
# Best fit parameters
boot.out31GA[which.max(boot.out31GA$Rn_max),]
# Conf.ints
quantile(boot.out31GA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out31GA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out31GA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
##################
# FULL BOOTSTRAP 2
#
# 10 runs take about 1 minute
for (i in 1:800)
{
# Sample and replace ALL monthly LFQ data
library(TropFishR)
library(parallel)
library(doParallel)
library(doRNG)
synLF5b <- synLFQ5
synLF5b$catch <- synLF5b$catch[,1:12]
synLF5b$dates <- synLF5b$dates[1:12]
synLF5b.t <- synLF5b
ncol(synLF5b.t$catch) # n = 12 montlhy LFQ data
mlen <- synLF5b.t$midLengths # midlengths
c.int <- (mlen[2]-mlen[1]) #class interval (1 cm)
# "m" ist the month(i.e, the catch file column), 1 to 25
for (m in 1:(ncol(synLF5b.t$catch))) {
catch.m <- synLF5b.t$catch[,m]
#plot(catch.m ~ mlen, xlim = c(20,90),main = "original LFD")# original LFD
#rebuild original data
reb.catch.m <- rep(mlen,catch.m) # recreates a proxy of original data
# take a sample
sam.reb.catch.m<- sample(reb.catch.m,replace = TRUE)
hist.sam.reb.catch.m <- hist(sam.reb.catch.m, breaks = c(mlen-(c.int/2)), plot = FALSE)
#replaces data in LFQ file ".t" with bootstrap sample data taken from ".t"
synLF5b.t$catch[1:length(hist.sam.reb.catch.m$counts) ,m] <- hist.sam.reb.catch.m$counts
}
synLF5b.RUN1 <- synLF5b.t # new sample
sample_seeds <- round(runif(100, max = 100, min = 1),0)
output2.new <- ELEFAN_GA(synLF5b.RUN1, seasonalised = F,
low_par = list(Linf = 40, K = 0.1, t_anchor = 0, C = 0, ts= 0),
up_par = list(Linf = 160, K = 7, t_anchor = 1, C = 0, ts = 1),
popSize = 50, maxiter = 20,
MA = 11, plot = F, plot.score = F , seed = sample(sample_seeds,1),
parallel = T, beep = F)
output2.new$par
output2.new$Rn_max
#plot(output2.new)
#date.time <- date()
Linf_out <- round(as.numeric(output2.new$par[1]), 2)
K_out <- round(as.numeric(output2.new$par[2]), 3)
t_anchor <- round(as.numeric(output2.new$par[3]), 3)
Phi_p <- round(as.numeric(output2.new$par[4]), 3)
Rn_max <- round(as.numeric(output2.new$Rn_max),3)
output.table <- data.frame(Linf_out, K_out, t_anchor, Phi_p,Rn_max )
write.table(output.table, file = "exp3BOOT_2_iter20_outp_table_GAdsyn5bh.csv",
sep = ";",col.names = FALSE, row.names = FALSE, append = TRUE )
rm(list=ls())
gc()
} # end of loop
boot.out32GA <- read.table( file = "exp3BOOT_2_iter20_outp_table_GAdsyn5bh.csv", sep = ";")
names(boot.out32GA) <- c("Linf", "K", "t_anchor", "Phi_p","Rn_max")
# Best fit parameters
boot.out32GA[which.max(boot.out32GA$Rn_max),]
# Conf.ints
quantile(boot.out32GA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out32GA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out32GA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
#######################################################
# Experiment 4
#
# Boot1 vs FULL Boot2, with _SA, Synlf5b
library(TropFishR)
# synLFQ5 # Linf = 80, K= 0.5,t_anchor = 0.25, C = 0, 2 years, perfect data
# use one year only
synLF5b <- synLFQ5
synLF5b$catch <- synLF5b$catch[,1:12]
synLF5b$dates <- synLF5b$dates[1:12]
##################
# BOOTSTRAP 1
# 10 runs take about 1 minute
# RUN BOOSTRAP 1 # maxit = 20, 10 runs = 1 MINUTE, sample_seeds <- round(runif(100, max = 100, min = 1),0)
for (i in 1:800) {
dat.time<- Sys.time() # Date and Time
dat.time.num<- as.numeric(dat.time) # Date and Time (numeric)
seed_glob <- i + dat.time.num # Global Seed value (date+time+"i"), used to Reset the Random Number Generator
set.seed(seed_glob) # Resets the Random Number Generator
library(TropFishR)
synLF5b <- synLFQ5
synLF5b$catch <- synLF5b$catch[,1:12]
synLF5b$dates <- synLF5b$dates[1:12]
Linf.out.vec <- rep(NA, 1)
K.out.vec <-rep(NA, 1)
t_anchor.out.vec <-rep(NA, 1)
phi.out.vec <- rep(NA, 1)
Rn_max.score.vec <- rep(NA, 1)
outp.table <- data.frame(Linf.out.vec, K.out.vec, t_anchor.out.vec, phi.out.vec, Rn_max.score.vec)
K.seeds_syn5b = round(runif(100, max = 3, min = 0.15),2)
Linf.seeds_syn5b = round(runif(100, max = 140, min = 50),1)
output <- ELEFAN_SA(synLF5b, seasonalised = F,
init_par = list(Linf = sample(Linf.seeds_syn5b,1), K = sample(K.seeds_syn5b,1), t_anchor = 0.5, C = 0, ts = 0),
low_par = list(Linf = 40, K = 0.1, t_anchor = 0, C = 0, ts= 0),
up_par = list(Linf = 160, K = 7, t_anchor = 1, C = 0, ts = 1),
beep = F, maxit = 20,
MA = 11, plot = F, plot.score = F
) # SA_time = 3
output$par
output$Rn_max
outp.table[1,1] <- round(as.numeric(output$par[1]), 2)
outp.table[1,2] <- round(as.numeric(output$par[2]), 3)
outp.table[1,3] <- round(as.numeric(output$par[3]), 3)
outp.table[1,4] <- round(as.numeric(output$par[4]), 3)
outp.table[1,5] <- round(output$Rn_max, 4)
#setwd("~/WG_New_Length-Based_Methods")
write.table(outp.table,file = "exp4BOOT_1_iter20_outp_table_SA_Syn5bh.csv" ,append = TRUE,
sep= ";", col.names = FALSE, row.names = FALSE)
# clear memory (avoids crashing the R session)
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
}
# setwd("~/WG_New_Length-Based_Methods")
boot.out41SA <- read.table(file = "exp4BOOT_1_iter20_outp_table_SA_Syn5bh.csv" , sep= ";")
names(boot.out41SA) <-c("Linf",
"K",
"t_anchor",
"Phi_p" ,
"Rn_max" )
# Best fit parameters
boot.out41SA[which.max(boot.out41SA$Rn_max),]
# Conf.ints
quantile(boot.out41SA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out41SA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out41SA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
##################
# FULL BOOTSTRAP 2
#
# 10 runs take about 1 minute
for (i in 1:800)
{
dat.time<- Sys.time() # Date and Time
dat.time.num<- as.numeric(dat.time) # Date and Time (numeric)
seed_glob <- i + dat.time.num # Global Seed value (date+time+"i"), used to Reset the Random Number Generator
set.seed(seed_glob) # Resets the Random Number Generator
# Sample and replace ALL monthly LFQ data
library(TropFishR)
synLF5b <- synLFQ5
synLF5b$catch <- synLF5b$catch[,1:12]
synLF5b$dates <- synLF5b$dates[1:12]
K.seeds_syn5b = round(runif(100, max = 3, min = 0.15),2)
Linf.seeds_syn5b = round(runif(100, max = 140, min = 50),1)
synLF5b.t <- synLF5b
ncol(synLF5b.t$catch) # n = x montlhy LFQ data
mlen <- synLF5b.t$midLengths # midlengths
c.int <- (mlen[2]-mlen[1]) #class interval (1 cm)
# "m" ist the month(i.e, the catch file column), 1 to 25
for (m in 1:(ncol(synLF5b.t$catch))) {
catch.m <- synLF5b.t$catch[,m]
#plot(catch.m ~ mlen, xlim = c(20,90),main = "original LFD")# original LFD
#rebuild original data
reb.catch.m <- rep(mlen,catch.m) # recreates a proxy of original data
# take a sample
sam.reb.catch.m<- sample(reb.catch.m,replace = TRUE)
hist.sam.reb.catch.m <- hist(sam.reb.catch.m, breaks = c(mlen-(c.int/2)), plot = FALSE)
#replaces data in LFQ file "2" with bootstrap sample data taken from "2"
synLF5b.t$catch[1:length(hist.sam.reb.catch.m$counts) ,m] <- hist.sam.reb.catch.m$counts
}
synLF5b.run1 <- synLF5b.t # new sample
output2.new <- ELEFAN_SA(synLF5b.run1, seasonalised = F,
init_par = list(Linf = sample(Linf.seeds_syn5b,1), K = sample(K.seeds_syn5b,1), t_anchor = 0.5, C = 0, ts = 0),
low_par = list(Linf = 40, K = 0.1, t_anchor = 0, C = 0, ts= 0),
up_par = list(Linf = 160, K = 7, t_anchor = 1, C = 0, ts = 1),
beep = F, maxit = 20,
MA = 11, plot = F, plot.score = F
)#SA_time = 3
output2.new$par
output2.new$Rn_max
output2.new$par
#plot(output2.new)
#date.time <- date()
Linf_out <- round(as.numeric(output2.new$par[1]), 2)
K_out <- round(as.numeric(output2.new$par[2]), 3)
t_anchor <- round(as.numeric(output2.new$par[3]), 3)
Phi_p <- round(as.numeric(output2.new$par[4]), 3)
Rn_max <- round(as.numeric(output2.new$Rn_max),3)
output.table <- data.frame(Linf_out, K_out, t_anchor, Phi_p,Rn_max )
write.table(output.table, file = "exp4BOOT_2_iter20_outp_table_SA_Syn5bh.csv", sep = ";",
col.names = FALSE, row.names = FALSE, append = TRUE )
rm(list=ls())
gc()
} # end of loop
# setwd("~/WG_New_Length-Based_Methods")
boot.out42SA <- read.table( file = "exp4BOOT_2_iter20_outp_table_SA_Syn5bh.csv", sep = ";")
names(boot.out42SA) <- c("Linf", "K", "t_anchor", "Phi_p","Rn_max")
# Best fit parameters
boot.out42SA[which.max(boot.out42SA$Rn_max),]
# Conf.ints
quantile(boot.out42SA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out42SA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out42SA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
#######################################################
# Experiment 5
#
# Boot1 vs FULL Boot2, with RSA_fun, alba
#
#function()
RSA_plot_for_LINUX_no_plot_b <- function (x, Linf_fix = NA, Linf_range = NA, K_range = exp(seq(log(0.1),
log(10), length.out = 100)), C = 0, ts = 0,MA = 9, addl.sqrt = FALSE,
agemax = NULL, flagging.out = TRUE, method = "optimise",
cross.date = NULL, cross.midLength = NULL, cross.max = FALSE,
hide.progressbar = FALSE, plot = FALSE, contour = FALSE,
add.values = TRUE, rsa.colors = terrain.colors(20), plot_title = TRUE)
{
res <- x
classes <- res$midLengths
catch <- res$catch
dates <- res$dates
interval <- (classes[2] - classes[1])/2
n_samples <- dim(catch)[2]
n_classes <- length(classes)
if (is.na(Linf_fix) & is.na(Linf_range[1]))
Linf_range <- seq(classes[n_classes] - 5, classes[n_classes] +
5, 1)
res <- lfqRestructure(res, MA = MA, addl.sqrt = addl.sqrt)
catch_aAF_F <- res$rcounts
peaks_mat <- res$peaks_mat
ASP <- res$ASP
if (!is.na(Linf_fix)) {
Linfs <- Linf_fix
}
else Linfs <- Linf_range
Ks <- K_range
sofun <- function(tanch, lfq, par, agemax, flagging.out) {
Lt <- lfqFitCurves(lfq, par = list(Linf = par[1], K = par[2],
t_anchor = tanch, C = par[4], ts = par[5]), flagging.out = flagging.out,
agemax = agemax)
return(Lt$ESP)
}
ESP_tanch_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
ESP_list_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
writeLines(paste("Optimisation procuedure of ELEFAN is running. \nThis will take some time. \nThe process bar will inform you about the process of the calculations.",
sep = " "))
flush.console()
if (method == "cross") {
if (cross.max) {
max.rcount <- which.max(res$rcounts)
COMB <- expand.grid(midLengths = rev(rev(res$midLengths)),
dates = res$dates)
cross.date <- COMB$dates[max.rcount]
cross.midLength <- COMB$midLengths[max.rcount]
}
else {
if (is.null(cross.date) | is.null(cross.midLength)) {
stop("Must define both 'cross.date' and 'cross.midLength' when 'cross.max' equals FALSE")
}
}
}
if (!hide.progressbar) {
nlk <- prod(dim(ESP_list_L))
pb <- txtProgressBar(min = 1, max = nlk, style = 3)
counter <- 1
}
for (li in 1:length(Linfs)) {
ESP_tanch_K <- rep(NA, length(Ks))
ESP_list_K <- rep(NA, length(Ks))
for (ki in 1:length(Ks)) {
if (is.null(agemax)) {
agemax.i <- ceiling((1/-Ks[ki]) * log(1 - ((Linfs[li] *
0.95)/Linfs[li])))
}
else {
agemax.i <- agemax
}
if (method == "cross") {
t0s <- seq(floor(min(date2yeardec(res$dates)) -
agemax.i), ceiling(max(date2yeardec(res$dates))),
by = 0.01)
Ltlut <- Linfs[li] * (1 - exp(-(Ks[ki] * (date2yeardec(cross.date) -
t0s) + (((C * Ks[ki])/(2 * pi)) * sin(2 *
pi * (date2yeardec(cross.date) - ts))) - (((C *
Ks[ki])/(2 * pi)) * sin(2 * pi * (t0s - ts))))))
t_anchor <- t0s[which.min((Ltlut - cross.midLength)^2)]%%1
resis <- sofun(tanch = t_anchor, lfq = res,
par = c(Linfs[li], Ks[ki], NA, C, ts), agemax = agemax.i,
flagging.out = flagging.out)
ESP_list_K[ki] <- resis
ESP_tanch_K[ki] <- t_anchor
}
if (method == "optimise") {
resis <- optimise(f = sofun, lower = 0, upper = 1,
lfq = res, par = c(Linfs[li], Ks[ki], NA,
C, ts), agemax = agemax.i, flagging.out = flagging.out,
tol = 0.001, maximum = TRUE)
ESP_list_K[ki] <- resis$objective
ESP_tanch_K[ki] <- resis$maximum
}
if (!hide.progressbar) {
setTxtProgressBar(pb, counter)
counter <- counter + 1
}
}
ESP_list_L[, li] <- ESP_list_K
ESP_tanch_L[, li] <- ESP_tanch_K
}
dimnames(ESP_list_L) <- list(Ks, Linfs)
score_mat <- round((10^(ESP_list_L/ASP))/10, digits = 3)
rownames(score_mat) <- round(as.numeric(rownames(score_mat)),
digits = 2)
dimnames(ESP_tanch_L) <- list(Ks, Linfs)
if (is.na(Linf_fix)) {
plot_dat <- reshape2::melt(score_mat)
#image(x = Linfs, y = Ks, z = t(score_mat), col = rsa.colors,
# ylab = "K", xlab = "Linf")
if (plot_title)
# title("Response surface analysis", line = 1)
if (contour) {
# contour(x = Linfs, y = Ks, z = t(score_mat), add = TRUE)
}
else {
if (is.numeric(contour)) {
#contour(x = Linfs, y = Ks, z = t(score_mat),
# add = TRUE, nlevels = contour)
}
else {
if (add.values) {
# text(x = plot_dat$Var2, y = plot_dat$Var1,
# round(as.numeric(plot_dat$value), digits = 2),
# cex = 0.6)
}
}
}
}
else {
if (all(Ks %in% exp(seq(log(0.1), log(10), length.out = 100)))) {
K_labels <- c(seq(0.1, 1, 0.1), 2:10)
K_plot <- log10(Ks)
K_ats <- log10(K_labels)
}
else {
K_labels <- Ks
K_plot <- Ks
K_ats <- Ks
}
phis <- round(log10(K_labels) + 2 * log10(Linfs), digits = 2)
op <- par(mar = c(12, 5, 4, 2))
# plot(K_plot, score_mat, type = "l", ylim = c(0, max(score_mat,
# na.rm = TRUE) * 1.4), ylab = "Score function", xlab = "Growth constant K (/year)",
# col = "red", lwd = 2, xaxt = "n")
#axis(1, at = K_ats, labels = K_labels)
#axis(1, at = K_ats, labels = phis, line = 5.5)
#mtext(text = expression(paste("Growth performance index (",
# phi, "')")), side = 1, line = 8.5)
#grid(nx = 0, NULL, lty = 6, col = "gray40")
#abline(v = K_ats, lty = 6, col = "gray40")
#if (plot_title)
# title("K-Scan", line = 1)
par(op)
}
Rn_max <- max(score_mat, na.rm = TRUE)[1]
idxs <- which(score_mat == Rn_max, arr.ind = TRUE)[1, ]
Linfest <- as.numeric(as.character(colnames(score_mat)[idxs[2]]))
Kest <- as.numeric(as.character(rownames(score_mat)[idxs[1]]))
tanchest <- as.numeric(as.character(ESP_tanch_L[idxs[1],
idxs[2]]))
phiL <- log10(Kest) + 2 * log10(Linfest)
pars <- list(Linf = Linfest, K = Kest, t_anchor = tanchest,
C = C, ts = ts, phiL = phiL)
final_res <- lfqFitCurves(lfq = res, par = pars, flagging.out = flagging.out,
agemax = agemax)
ret <- c(res, list(score_mat = score_mat, t_anchor_mat = ESP_tanch_L,
ncohort = final_res$ncohort, agemax = final_res$agemax,
par = pars, Rn_max = Rn_max))
class(ret) <- "lfq"
if (plot) {
#plot(ret, Fname = "rcounts")
Lt <- lfqFitCurves(ret, par = pars, draw = TRUE)
}
return(ret)
}
#
##################
# BOOTSTRAP 1
# 10 runs take about 1 minute
test.B1runs <- 1
# Boot 1 with RSA_plot fit should produce N runs with identical estimates (same data -> same result)
# RUN BOOSTRAP 1 # 5 runs = 1/2 MINUTE , no seeds, no maxit, but: effect of MA?
for (i in 1:test.B1runs) {
MA.input = 9
library(TropFishR)
Linf.out.vec <- rep(NA, 1)
K.out.vec <-rep(NA, 1)
t_anchor.out.vec <-rep(NA, 1)
phi.out.vec <- rep(NA, 1)
Rn_max.score.vec <- rep(NA, 1)
outp.table <- data.frame(Linf.out.vec, K.out.vec, t_anchor.out.vec, phi.out.vec, Rn_max.score.vec,
MA.input)
alba2 <- readRDS("alba2.rds")
RSA_plot_for_LINUX_no_plot_b <- function (x, Linf_fix = NA, Linf_range = NA, K_range = exp(seq(log(0.1),
log(10), length.out = 100)), C = 0, ts = 0,MA = 9, addl.sqrt = FALSE,
agemax = NULL, flagging.out = TRUE, method = "optimise",
cross.date = NULL, cross.midLength = NULL, cross.max = FALSE,
hide.progressbar = FALSE, plot = FALSE, contour = FALSE,
add.values = TRUE, rsa.colors = terrain.colors(20), plot_title = TRUE)
{
res <- x
classes <- res$midLengths
catch <- res$catch
dates <- res$dates
interval <- (classes[2] - classes[1])/2
n_samples <- dim(catch)[2]
n_classes <- length(classes)
if (is.na(Linf_fix) & is.na(Linf_range[1]))
Linf_range <- seq(classes[n_classes] - 5, classes[n_classes] +
5, 1)
res <- lfqRestructure(res, MA = MA, addl.sqrt = addl.sqrt)
catch_aAF_F <- res$rcounts
peaks_mat <- res$peaks_mat
ASP <- res$ASP
if (!is.na(Linf_fix)) {
Linfs <- Linf_fix
}
else Linfs <- Linf_range
Ks <- K_range
sofun <- function(tanch, lfq, par, agemax, flagging.out) {
Lt <- lfqFitCurves(lfq, par = list(Linf = par[1], K = par[2],
t_anchor = tanch, C = par[4], ts = par[5]), flagging.out = flagging.out,
agemax = agemax)
return(Lt$ESP)
}
ESP_tanch_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
ESP_list_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
writeLines(paste("Optimisation procuedure of ELEFAN is running. \nThis will take some time. \nThe process bar will inform you about the process of the calculations.",
sep = " "))
flush.console()
if (method == "cross") {
if (cross.max) {
max.rcount <- which.max(res$rcounts)
COMB <- expand.grid(midLengths = rev(rev(res$midLengths)),
dates = res$dates)
cross.date <- COMB$dates[max.rcount]
cross.midLength <- COMB$midLengths[max.rcount]
}
else {
if (is.null(cross.date) | is.null(cross.midLength)) {
stop("Must define both 'cross.date' and 'cross.midLength' when 'cross.max' equals FALSE")
}
}
}
if (!hide.progressbar) {
nlk <- prod(dim(ESP_list_L))
pb <- txtProgressBar(min = 1, max = nlk, style = 3)
counter <- 1
}
for (li in 1:length(Linfs)) {
ESP_tanch_K <- rep(NA, length(Ks))
ESP_list_K <- rep(NA, length(Ks))
for (ki in 1:length(Ks)) {
if (is.null(agemax)) {
agemax.i <- ceiling((1/-Ks[ki]) * log(1 - ((Linfs[li] *
0.95)/Linfs[li])))
}
else {
agemax.i <- agemax
}
if (method == "cross") {
t0s <- seq(floor(min(date2yeardec(res$dates)) -
agemax.i), ceiling(max(date2yeardec(res$dates))),
by = 0.01)
Ltlut <- Linfs[li] * (1 - exp(-(Ks[ki] * (date2yeardec(cross.date) -
t0s) + (((C * Ks[ki])/(2 * pi)) * sin(2 *
pi * (date2yeardec(cross.date) - ts))) - (((C *
Ks[ki])/(2 * pi)) * sin(2 * pi * (t0s - ts))))))
t_anchor <- t0s[which.min((Ltlut - cross.midLength)^2)]%%1
resis <- sofun(tanch = t_anchor, lfq = res,
par = c(Linfs[li], Ks[ki], NA, C, ts), agemax = agemax.i,
flagging.out = flagging.out)
ESP_list_K[ki] <- resis
ESP_tanch_K[ki] <- t_anchor
}
if (method == "optimise") {
resis <- optimise(f = sofun, lower = 0, upper = 1,
lfq = res, par = c(Linfs[li], Ks[ki], NA,
C, ts), agemax = agemax.i, flagging.out = flagging.out,
tol = 0.001, maximum = TRUE)
ESP_list_K[ki] <- resis$objective
ESP_tanch_K[ki] <- resis$maximum
}
if (!hide.progressbar) {
setTxtProgressBar(pb, counter)
counter <- counter + 1
}
}
ESP_list_L[, li] <- ESP_list_K
ESP_tanch_L[, li] <- ESP_tanch_K
}
dimnames(ESP_list_L) <- list(Ks, Linfs)
score_mat <- round((10^(ESP_list_L/ASP))/10, digits = 3)
rownames(score_mat) <- round(as.numeric(rownames(score_mat)),
digits = 2)
dimnames(ESP_tanch_L) <- list(Ks, Linfs)
if (is.na(Linf_fix)) {
plot_dat <- reshape2::melt(score_mat)
#image(x = Linfs, y = Ks, z = t(score_mat), col = rsa.colors,
# ylab = "K", xlab = "Linf")
if (plot_title)
# title("Response surface analysis", line = 1)
if (contour) {
# contour(x = Linfs, y = Ks, z = t(score_mat), add = TRUE)
}
else {
if (is.numeric(contour)) {
#contour(x = Linfs, y = Ks, z = t(score_mat),
# add = TRUE, nlevels = contour)
}
else {
if (add.values) {
# text(x = plot_dat$Var2, y = plot_dat$Var1,
# round(as.numeric(plot_dat$value), digits = 2),
# cex = 0.6)
}
}
}
}
else {
if (all(Ks %in% exp(seq(log(0.1), log(10), length.out = 100)))) {
K_labels <- c(seq(0.1, 1, 0.1), 2:10)
K_plot <- log10(Ks)
K_ats <- log10(K_labels)
}
else {
K_labels <- Ks
K_plot <- Ks
K_ats <- Ks
}
phis <- round(log10(K_labels) + 2 * log10(Linfs), digits = 2)
op <- par(mar = c(12, 5, 4, 2))
# plot(K_plot, score_mat, type = "l", ylim = c(0, max(score_mat,
# na.rm = TRUE) * 1.4), ylab = "Score function", xlab = "Growth constant K (/year)",
# col = "red", lwd = 2, xaxt = "n")
#axis(1, at = K_ats, labels = K_labels)
#axis(1, at = K_ats, labels = phis, line = 5.5)
#mtext(text = expression(paste("Growth performance index (",
# phi, "')")), side = 1, line = 8.5)
#grid(nx = 0, NULL, lty = 6, col = "gray40")
#abline(v = K_ats, lty = 6, col = "gray40")
#if (plot_title)
# title("K-Scan", line = 1)
par(op)
}
Rn_max <- max(score_mat, na.rm = TRUE)[1]
idxs <- which(score_mat == Rn_max, arr.ind = TRUE)[1, ]
Linfest <- as.numeric(as.character(colnames(score_mat)[idxs[2]]))
Kest <- as.numeric(as.character(rownames(score_mat)[idxs[1]]))
tanchest <- as.numeric(as.character(ESP_tanch_L[idxs[1],
idxs[2]]))
phiL <- log10(Kest) + 2 * log10(Linfest)
pars <- list(Linf = Linfest, K = Kest, t_anchor = tanchest,
C = C, ts = ts, phiL = phiL)
final_res <- lfqFitCurves(lfq = res, par = pars, flagging.out = flagging.out,
agemax = agemax)
ret <- c(res, list(score_mat = score_mat, t_anchor_mat = ESP_tanch_L,
ncohort = final_res$ncohort, agemax = final_res$agemax,
par = pars, Rn_max = Rn_max))
class(ret) <- "lfq"
if (plot) {
#plot(ret, Fname = "rcounts")
Lt <- lfqFitCurves(ret, par = pars, draw = TRUE)
}
return(ret)
}
#
output <- RSA_plot_for_LINUX_no_plot_b(alba2, Linf_range = seq(8,13,0.5),
K_range = seq(0.1,3,0.2),
C = 0, ts = 0, MA = MA.input,
contour = 3, plot = FALSE )
output$par
output$Rn_max
outp.table[1,1] <- round(as.numeric(output$par[1]), 2)
outp.table[1,2] <- round(as.numeric(output$par[2]), 3)
outp.table[1,3] <- round(as.numeric(output$par[3]), 3)
outp.table[1,4] <- round(as.numeric(output$par[6]), 3)
outp.table[1,5] <- round(output$Rn_max, 4)
#setwd("~/WG_New_Length-Based_Methods")
write.table(outp.table,file = "exp5BOOT_1_outp_table_RSAcalbaMA9.csv" ,
append = TRUE,sep= ";", col.names = FALSE, row.names = FALSE)
# clear memory (avoids crashing the R session)
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
}
# setwd("~/WG_New_Length-Based_Methods")
boot.out51RSA <- read.table(file = "exp5BOOT_1_outp_table_RSAcalbaMA9.csv" , sep= ";")
names(boot.out51RSA) <-c("Linf",
"K",
"t_anchor",
"Phi_p" ,
"Rn_max" ,
"MA")
# Best fit parameters
boot.out11GA[which.max(boot.out11GA$Rn_max),]
# Conf.ints
quantile(boot.out11GA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out11GA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out11GA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
################# RSA plot (1 plot, not for BOOTSTRAP)
# For one-time use only
# Big picture and extremley slow RSA (wide range, high resolution) - not for BOOTSTRAP - just for one plot
output.RSA.mat.for.plotMA5 <- RSA_plot_for_LINUX_no_plot_b(alba2, Linf_range = seq(6,18,0.1),
K_range = seq(0.4,7,0.1),
C = 0, ts = 0,MA = 9,
contour = 3, plot = F )
# "best fit" Rn values: par$Linf = 10.3, par$K = 1.4
saveRDS(output.RSA.mat.for.plot,file = "output_RSA_mat_for_plotMA9.rds")
rm(output.RSA.mat.for.plot)
output.RSA.mat.for.plot <- readRDS("output_RSA_mat_for_plotMA9.rds")
output.RSA.mat.for.plot
# View(as.data.frame(output.RSA.mat.for.plot$score_mat))
# plot(output.RSA.mat.for.plot$score_mat)
##################
# FULL BOOTSTRAP 2 (narrow range, low resolution, no plot "fast" RSA plot)
#
# 10 runs take about 1 minute
for (i in 1:3)
{
# Sample and replace ALL monthly LFQ data
library(TropFishR)
alba2 <- readRDS("alba2.rds")
alba2.t <- alba2
ncol(alba2.t$catch) # n = x montlhy LFQ data
mlen <- alba2.t$midLengths # midlengths
c.int <- (mlen[2]-mlen[1]) #class interval (1 cm)
# "m" ist the month(i.e, the catch file column), 1 to 25
for (m in 1:(ncol(alba2.t$catch))) {
catch.m <- alba2.t$catch[,m]
#plot(catch.m ~ mlen, xlim = c(20,90),main = "original LFD")# original LFD
#rebuild original data
reb.catch.m <- rep(mlen,catch.m) # recreates a proxy of original data
# take a sample
sam.reb.catch.m<- sample(reb.catch.m,replace = TRUE)
hist.sam.reb.catch.m <- hist(sam.reb.catch.m, breaks = c(mlen-(c.int/2)), plot = FALSE)
#replaces data in LFQ file "2" with bootstrap sample data taken from "2"
alba2.t$catch[1:length(hist.sam.reb.catch.m$counts) ,m] <- hist.sam.reb.catch.m$counts
}
alba2.run1 <- alba2.t # new sample
RSA_plot_for_LINUX_no_plot_b <- function (x, Linf_fix = NA, Linf_range = NA, K_range = exp(seq(log(0.1),
log(10), length.out = 100)), C = 0, ts = 0,MA = 9, addl.sqrt = FALSE,
agemax = NULL, flagging.out = TRUE, method = "optimise",
cross.date = NULL, cross.midLength = NULL, cross.max = FALSE,
hide.progressbar = FALSE, plot = FALSE, contour = FALSE,
add.values = TRUE, rsa.colors = terrain.colors(20), plot_title = TRUE)
{
res <- x
classes <- res$midLengths
catch <- res$catch
dates <- res$dates
interval <- (classes[2] - classes[1])/2
n_samples <- dim(catch)[2]
n_classes <- length(classes)
if (is.na(Linf_fix) & is.na(Linf_range[1]))
Linf_range <- seq(classes[n_classes] - 5, classes[n_classes] +
5, 1)
res <- lfqRestructure(res, MA = MA, addl.sqrt = addl.sqrt)
catch_aAF_F <- res$rcounts
peaks_mat <- res$peaks_mat
ASP <- res$ASP
if (!is.na(Linf_fix)) {
Linfs <- Linf_fix
}
else Linfs <- Linf_range
Ks <- K_range
sofun <- function(tanch, lfq, par, agemax, flagging.out) {
Lt <- lfqFitCurves(lfq, par = list(Linf = par[1], K = par[2],
t_anchor = tanch, C = par[4], ts = par[5]), flagging.out = flagging.out,
agemax = agemax)
return(Lt$ESP)
}
ESP_tanch_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
ESP_list_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
writeLines(paste("Optimisation procuedure of ELEFAN is running. \nThis will take some time. \nThe process bar will inform you about the process of the calculations.",
sep = " "))
flush.console()
if (method == "cross") {
if (cross.max) {
max.rcount <- which.max(res$rcounts)
COMB <- expand.grid(midLengths = rev(rev(res$midLengths)),
dates = res$dates)
cross.date <- COMB$dates[max.rcount]
cross.midLength <- COMB$midLengths[max.rcount]
}
else {
if (is.null(cross.date) | is.null(cross.midLength)) {
stop("Must define both 'cross.date' and 'cross.midLength' when 'cross.max' equals FALSE")
}
}
}
if (!hide.progressbar) {
nlk <- prod(dim(ESP_list_L))
pb <- txtProgressBar(min = 1, max = nlk, style = 3)
counter <- 1
}
for (li in 1:length(Linfs)) {
ESP_tanch_K <- rep(NA, length(Ks))
ESP_list_K <- rep(NA, length(Ks))
for (ki in 1:length(Ks)) {
if (is.null(agemax)) {
agemax.i <- ceiling((1/-Ks[ki]) * log(1 - ((Linfs[li] *
0.95)/Linfs[li])))
}
else {
agemax.i <- agemax
}
if (method == "cross") {
t0s <- seq(floor(min(date2yeardec(res$dates)) -
agemax.i), ceiling(max(date2yeardec(res$dates))),
by = 0.01)
Ltlut <- Linfs[li] * (1 - exp(-(Ks[ki] * (date2yeardec(cross.date) -
t0s) + (((C * Ks[ki])/(2 * pi)) * sin(2 *
pi * (date2yeardec(cross.date) - ts))) - (((C *
Ks[ki])/(2 * pi)) * sin(2 * pi * (t0s - ts))))))
t_anchor <- t0s[which.min((Ltlut - cross.midLength)^2)]%%1
resis <- sofun(tanch = t_anchor, lfq = res,
par = c(Linfs[li], Ks[ki], NA, C, ts), agemax = agemax.i,
flagging.out = flagging.out)
ESP_list_K[ki] <- resis
ESP_tanch_K[ki] <- t_anchor
}
if (method == "optimise") {
resis <- optimise(f = sofun, lower = 0, upper = 1,
lfq = res, par = c(Linfs[li], Ks[ki], NA,
C, ts), agemax = agemax.i, flagging.out = flagging.out,
tol = 0.001, maximum = TRUE)
ESP_list_K[ki] <- resis$objective
ESP_tanch_K[ki] <- resis$maximum
}
if (!hide.progressbar) {
setTxtProgressBar(pb, counter)
counter <- counter + 1
}
}
ESP_list_L[, li] <- ESP_list_K
ESP_tanch_L[, li] <- ESP_tanch_K
}
dimnames(ESP_list_L) <- list(Ks, Linfs)
score_mat <- round((10^(ESP_list_L/ASP))/10, digits = 3)
rownames(score_mat) <- round(as.numeric(rownames(score_mat)),
digits = 2)
dimnames(ESP_tanch_L) <- list(Ks, Linfs)
if (is.na(Linf_fix)) {
plot_dat <- reshape2::melt(score_mat)
#image(x = Linfs, y = Ks, z = t(score_mat), col = rsa.colors,
# ylab = "K", xlab = "Linf")
if (plot_title)
# title("Response surface analysis", line = 1)
if (contour) {
# contour(x = Linfs, y = Ks, z = t(score_mat), add = TRUE)
}
else {
if (is.numeric(contour)) {
#contour(x = Linfs, y = Ks, z = t(score_mat),
# add = TRUE, nlevels = contour)
}
else {
if (add.values) {
# text(x = plot_dat$Var2, y = plot_dat$Var1,
# round(as.numeric(plot_dat$value), digits = 2),
# cex = 0.6)
}
}
}
}
else {
if (all(Ks %in% exp(seq(log(0.1), log(10), length.out = 100)))) {
K_labels <- c(seq(0.1, 1, 0.1), 2:10)
K_plot <- log10(Ks)
K_ats <- log10(K_labels)
}
else {
K_labels <- Ks
K_plot <- Ks
K_ats <- Ks
}
phis <- round(log10(K_labels) + 2 * log10(Linfs), digits = 2)
op <- par(mar = c(12, 5, 4, 2))
# plot(K_plot, score_mat, type = "l", ylim = c(0, max(score_mat,
# na.rm = TRUE) * 1.4), ylab = "Score function", xlab = "Growth constant K (/year)",
# col = "red", lwd = 2, xaxt = "n")
#axis(1, at = K_ats, labels = K_labels)
#axis(1, at = K_ats, labels = phis, line = 5.5)
#mtext(text = expression(paste("Growth performance index (",
# phi, "')")), side = 1, line = 8.5)
#grid(nx = 0, NULL, lty = 6, col = "gray40")
#abline(v = K_ats, lty = 6, col = "gray40")
#if (plot_title)
# title("K-Scan", line = 1)
par(op)
}
Rn_max <- max(score_mat, na.rm = TRUE)[1]
idxs <- which(score_mat == Rn_max, arr.ind = TRUE)[1, ]
Linfest <- as.numeric(as.character(colnames(score_mat)[idxs[2]]))
Kest <- as.numeric(as.character(rownames(score_mat)[idxs[1]]))
tanchest <- as.numeric(as.character(ESP_tanch_L[idxs[1],
idxs[2]]))
phiL <- log10(Kest) + 2 * log10(Linfest)
pars <- list(Linf = Linfest, K = Kest, t_anchor = tanchest,
C = C, ts = ts, phiL = phiL)
final_res <- lfqFitCurves(lfq = res, par = pars, flagging.out = flagging.out,
agemax = agemax)
ret <- c(res, list(score_mat = score_mat, t_anchor_mat = ESP_tanch_L,
ncohort = final_res$ncohort, agemax = final_res$agemax,
par = pars, Rn_max = Rn_max))
class(ret) <- "lfq"
if (plot) {
#plot(ret, Fname = "rcounts")
Lt <- lfqFitCurves(ret, par = pars, draw = TRUE)
}
return(ret)
}
# Big picture and veeery slow RSA (wide range, high resolution)
#output2.new <- RSA_plot_for_LINUX_no_plot_b(alba2.run1, Linf_range = seq(7,18,0.2),
# K_range = seq(0.4,7,0.2),
# C = 0, ts = 0, MA = 9,
# contour = 3, plot = FALSE )
# narrow and fast RSA (searches whithin narrow range)
output2.new <- RSA_plot_for_LINUX_no_plot_b(alba2.run1, Linf_range = seq(8,13,0.2),
K_range = seq(1,7,0.2),
C = 0, ts = 0, MA = 9,
contour = 3, plot = FALSE )
output2.new$par
output2.new$Rn_max
output2.new$par
#plot(output2.new)
#date.time <- date()
Linf_out <- round(as.numeric(output2.new$par[1]), 2)
K_out <- round(as.numeric(output2.new$par[2]), 3)
t_anchor <- round(as.numeric(output2.new$par[3]), 3)
Phi_p <- round(as.numeric(output2.new$par[6]), 3)
Rn_max <- round(as.numeric(output2.new$Rn_max),3)
output.table <- data.frame(Linf_out, K_out, t_anchor, Phi_p,Rn_max)
write.table(output.table, file = "exp5BOOT_2__outp_table_RSAalbaMA9.csv",
sep = ";",col.names = FALSE, row.names = FALSE, append = TRUE )
rm(list=ls())
gc()
} # end of loop
boot.out52RSA_alba <- read.table( file = "exp5BOOT_2__outp_table_RSAalbaMA9.csv", sep = ";")
names(boot.out52RSA_alba) <- c("Linf", "K", "t_anchor", "Phi_p","Rn_max")
boot.out52RSA_alba
# Best fit parameters
boot.out52RSA_alba[which.max(boot.out52RSA_alba$Rn_max),]
# Conf.ints
quantile(boot.out52RSA_alba$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out52RSA_alba$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out52RSA_alba$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
##
#######################################################
# Experiment 6
#
# Boot1 vs FULL Boot2, with RSA_fun, synlfq5b
#
#function()
RSA_plot_for_LINUX_no_plot_b <- function (x, Linf_fix = NA, Linf_range = NA, K_range = exp(seq(log(0.1),
log(10), length.out = 100)), C = 0, ts = 0,MA = 9, addl.sqrt = FALSE,
agemax = NULL, flagging.out = TRUE, method = "optimise",
cross.date = NULL, cross.midLength = NULL, cross.max = FALSE,
hide.progressbar = FALSE, plot = FALSE, contour = FALSE,
add.values = TRUE, rsa.colors = terrain.colors(20), plot_title = TRUE)
{
res <- x
classes <- res$midLengths
catch <- res$catch
dates <- res$dates
interval <- (classes[2] - classes[1])/2
n_samples <- dim(catch)[2]
n_classes <- length(classes)
if (is.na(Linf_fix) & is.na(Linf_range[1]))
Linf_range <- seq(classes[n_classes] - 5, classes[n_classes] +
5, 1)
res <- lfqRestructure(res, MA = MA, addl.sqrt = addl.sqrt)
catch_aAF_F <- res$rcounts
peaks_mat <- res$peaks_mat
ASP <- res$ASP
if (!is.na(Linf_fix)) {
Linfs <- Linf_fix
}
else Linfs <- Linf_range
Ks <- K_range
sofun <- function(tanch, lfq, par, agemax, flagging.out) {
Lt <- lfqFitCurves(lfq, par = list(Linf = par[1], K = par[2],
t_anchor = tanch, C = par[4], ts = par[5]), flagging.out = flagging.out,
agemax = agemax)
return(Lt$ESP)
}
ESP_tanch_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
ESP_list_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
writeLines(paste("Optimisation procuedure of ELEFAN is running. \nThis will take some time. \nThe process bar will inform you about the process of the calculations.",
sep = " "))
flush.console()
if (method == "cross") {
if (cross.max) {
max.rcount <- which.max(res$rcounts)
COMB <- expand.grid(midLengths = rev(rev(res$midLengths)),
dates = res$dates)
cross.date <- COMB$dates[max.rcount]
cross.midLength <- COMB$midLengths[max.rcount]
}
else {
if (is.null(cross.date) | is.null(cross.midLength)) {
stop("Must define both 'cross.date' and 'cross.midLength' when 'cross.max' equals FALSE")
}
}
}
if (!hide.progressbar) {
nlk <- prod(dim(ESP_list_L))
pb <- txtProgressBar(min = 1, max = nlk, style = 3)
counter <- 1
}
for (li in 1:length(Linfs)) {
ESP_tanch_K <- rep(NA, length(Ks))
ESP_list_K <- rep(NA, length(Ks))
for (ki in 1:length(Ks)) {
if (is.null(agemax)) {
agemax.i <- ceiling((1/-Ks[ki]) * log(1 - ((Linfs[li] *
0.95)/Linfs[li])))
}
else {
agemax.i <- agemax
}
if (method == "cross") {
t0s <- seq(floor(min(date2yeardec(res$dates)) -
agemax.i), ceiling(max(date2yeardec(res$dates))),
by = 0.01)
Ltlut <- Linfs[li] * (1 - exp(-(Ks[ki] * (date2yeardec(cross.date) -
t0s) + (((C * Ks[ki])/(2 * pi)) * sin(2 *
pi * (date2yeardec(cross.date) - ts))) - (((C *
Ks[ki])/(2 * pi)) * sin(2 * pi * (t0s - ts))))))
t_anchor <- t0s[which.min((Ltlut - cross.midLength)^2)]%%1
resis <- sofun(tanch = t_anchor, lfq = res,
par = c(Linfs[li], Ks[ki], NA, C, ts), agemax = agemax.i,
flagging.out = flagging.out)
ESP_list_K[ki] <- resis
ESP_tanch_K[ki] <- t_anchor
}
if (method == "optimise") {
resis <- optimise(f = sofun, lower = 0, upper = 1,
lfq = res, par = c(Linfs[li], Ks[ki], NA,
C, ts), agemax = agemax.i, flagging.out = flagging.out,
tol = 0.001, maximum = TRUE)
ESP_list_K[ki] <- resis$objective
ESP_tanch_K[ki] <- resis$maximum
}
if (!hide.progressbar) {
setTxtProgressBar(pb, counter)
counter <- counter + 1
}
}
ESP_list_L[, li] <- ESP_list_K
ESP_tanch_L[, li] <- ESP_tanch_K
}
dimnames(ESP_list_L) <- list(Ks, Linfs)
score_mat <- round((10^(ESP_list_L/ASP))/10, digits = 3)
rownames(score_mat) <- round(as.numeric(rownames(score_mat)),
digits = 2)
dimnames(ESP_tanch_L) <- list(Ks, Linfs)
if (is.na(Linf_fix)) {
plot_dat <- reshape2::melt(score_mat)
#image(x = Linfs, y = Ks, z = t(score_mat), col = rsa.colors,
# ylab = "K", xlab = "Linf")
if (plot_title)
# title("Response surface analysis", line = 1)
if (contour) {
# contour(x = Linfs, y = Ks, z = t(score_mat), add = TRUE)
}
else {
if (is.numeric(contour)) {
#contour(x = Linfs, y = Ks, z = t(score_mat),
# add = TRUE, nlevels = contour)
}
else {
if (add.values) {
# text(x = plot_dat$Var2, y = plot_dat$Var1,
# round(as.numeric(plot_dat$value), digits = 2),
# cex = 0.6)
}
}
}
}
else {
if (all(Ks %in% exp(seq(log(0.1), log(10), length.out = 100)))) {
K_labels <- c(seq(0.1, 1, 0.1), 2:10)
K_plot <- log10(Ks)
K_ats <- log10(K_labels)
}
else {
K_labels <- Ks
K_plot <- Ks
K_ats <- Ks
}
phis <- round(log10(K_labels) + 2 * log10(Linfs), digits = 2)
op <- par(mar = c(12, 5, 4, 2))
# plot(K_plot, score_mat, type = "l", ylim = c(0, max(score_mat,
# na.rm = TRUE) * 1.4), ylab = "Score function", xlab = "Growth constant K (/year)",
# col = "red", lwd = 2, xaxt = "n")
#axis(1, at = K_ats, labels = K_labels)
#axis(1, at = K_ats, labels = phis, line = 5.5)
#mtext(text = expression(paste("Growth performance index (",
# phi, "')")), side = 1, line = 8.5)
#grid(nx = 0, NULL, lty = 6, col = "gray40")
#abline(v = K_ats, lty = 6, col = "gray40")
#if (plot_title)
# title("K-Scan", line = 1)
par(op)
}
Rn_max <- max(score_mat, na.rm = TRUE)[1]
idxs <- which(score_mat == Rn_max, arr.ind = TRUE)[1, ]
Linfest <- as.numeric(as.character(colnames(score_mat)[idxs[2]]))
Kest <- as.numeric(as.character(rownames(score_mat)[idxs[1]]))
tanchest <- as.numeric(as.character(ESP_tanch_L[idxs[1],
idxs[2]]))
phiL <- log10(Kest) + 2 * log10(Linfest)
pars <- list(Linf = Linfest, K = Kest, t_anchor = tanchest,
C = C, ts = ts, phiL = phiL)
final_res <- lfqFitCurves(lfq = res, par = pars, flagging.out = flagging.out,
agemax = agemax)
ret <- c(res, list(score_mat = score_mat, t_anchor_mat = ESP_tanch_L,
ncohort = final_res$ncohort, agemax = final_res$agemax,
par = pars, Rn_max = Rn_max))
class(ret) <- "lfq"
if (plot) {
#plot(ret, Fname = "rcounts")
Lt <- lfqFitCurves(ret, par = pars, draw = TRUE)
}
return(ret)
}
#
##################
# BOOTSTRAP 1
# 10 runs take about 1 minute
test.B1runs <- 1
# Boot 1 with RSA_plot fit should produce N runs with identical estimates (same data -> same result)
# RUN BOOSTRAP 1 # 5 runs = 1/2 MINUTE , no seeds, no maxit, but: effect of MA?
for (i in 1:test.B1runs) {
MA.input = 9
library(TropFishR)
# use one year only
synLF5b <- synLFQ5
synLF5b$catch <- synLF5b$catch[,1:12]
synLF5b$dates <- synLF5b$dates[1:12]
Linf.out.vec <- rep(NA, 1)
K.out.vec <-rep(NA, 1)
t_anchor.out.vec <-rep(NA, 1)
phi.out.vec <- rep(NA, 1)
Rn_max.score.vec <- rep(NA, 1)
outp.table <- data.frame(Linf.out.vec, K.out.vec, t_anchor.out.vec, phi.out.vec, Rn_max.score.vec,
MA.input)
RSA_plot_for_LINUX_no_plot_b <- function (x, Linf_fix = NA, Linf_range = NA, K_range = exp(seq(log(0.1),
log(10), length.out = 100)), C = 0, ts = 0,MA = 9, addl.sqrt = FALSE,
agemax = NULL, flagging.out = TRUE, method = "optimise",
cross.date = NULL, cross.midLength = NULL, cross.max = FALSE,
hide.progressbar = FALSE, plot = FALSE, contour = FALSE,
add.values = TRUE, rsa.colors = terrain.colors(20), plot_title = TRUE)
{
res <- x
classes <- res$midLengths
catch <- res$catch
dates <- res$dates
interval <- (classes[2] - classes[1])/2
n_samples <- dim(catch)[2]
n_classes <- length(classes)
if (is.na(Linf_fix) & is.na(Linf_range[1]))
Linf_range <- seq(classes[n_classes] - 5, classes[n_classes] +
5, 1)
res <- lfqRestructure(res, MA = MA, addl.sqrt = addl.sqrt)
catch_aAF_F <- res$rcounts
peaks_mat <- res$peaks_mat
ASP <- res$ASP
if (!is.na(Linf_fix)) {
Linfs <- Linf_fix
}
else Linfs <- Linf_range
Ks <- K_range
sofun <- function(tanch, lfq, par, agemax, flagging.out) {
Lt <- lfqFitCurves(lfq, par = list(Linf = par[1], K = par[2],
t_anchor = tanch, C = par[4], ts = par[5]), flagging.out = flagging.out,
agemax = agemax)
return(Lt$ESP)
}
ESP_tanch_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
ESP_list_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
writeLines(paste("Optimisation procuedure of ELEFAN is running. \nThis will take some time. \nThe process bar will inform you about the process of the calculations.",
sep = " "))
flush.console()
if (method == "cross") {
if (cross.max) {
max.rcount <- which.max(res$rcounts)
COMB <- expand.grid(midLengths = rev(rev(res$midLengths)),
dates = res$dates)
cross.date <- COMB$dates[max.rcount]
cross.midLength <- COMB$midLengths[max.rcount]
}
else {
if (is.null(cross.date) | is.null(cross.midLength)) {
stop("Must define both 'cross.date' and 'cross.midLength' when 'cross.max' equals FALSE")
}
}
}
if (!hide.progressbar) {
nlk <- prod(dim(ESP_list_L))
pb <- txtProgressBar(min = 1, max = nlk, style = 3)
counter <- 1
}
for (li in 1:length(Linfs)) {
ESP_tanch_K <- rep(NA, length(Ks))
ESP_list_K <- rep(NA, length(Ks))
for (ki in 1:length(Ks)) {
if (is.null(agemax)) {
agemax.i <- ceiling((1/-Ks[ki]) * log(1 - ((Linfs[li] *
0.95)/Linfs[li])))
}
else {
agemax.i <- agemax
}
if (method == "cross") {
t0s <- seq(floor(min(date2yeardec(res$dates)) -
agemax.i), ceiling(max(date2yeardec(res$dates))),
by = 0.01)
Ltlut <- Linfs[li] * (1 - exp(-(Ks[ki] * (date2yeardec(cross.date) -
t0s) + (((C * Ks[ki])/(2 * pi)) * sin(2 *
pi * (date2yeardec(cross.date) - ts))) - (((C *
Ks[ki])/(2 * pi)) * sin(2 * pi * (t0s - ts))))))
t_anchor <- t0s[which.min((Ltlut - cross.midLength)^2)]%%1
resis <- sofun(tanch = t_anchor, lfq = res,
par = c(Linfs[li], Ks[ki], NA, C, ts), agemax = agemax.i,
flagging.out = flagging.out)
ESP_list_K[ki] <- resis
ESP_tanch_K[ki] <- t_anchor
}
if (method == "optimise") {
resis <- optimise(f = sofun, lower = 0, upper = 1,
lfq = res, par = c(Linfs[li], Ks[ki], NA,
C, ts), agemax = agemax.i, flagging.out = flagging.out,
tol = 0.001, maximum = TRUE)
ESP_list_K[ki] <- resis$objective
ESP_tanch_K[ki] <- resis$maximum
}
if (!hide.progressbar) {
setTxtProgressBar(pb, counter)
counter <- counter + 1
}
}
ESP_list_L[, li] <- ESP_list_K
ESP_tanch_L[, li] <- ESP_tanch_K
}
dimnames(ESP_list_L) <- list(Ks, Linfs)
score_mat <- round((10^(ESP_list_L/ASP))/10, digits = 3)
rownames(score_mat) <- round(as.numeric(rownames(score_mat)),
digits = 2)
dimnames(ESP_tanch_L) <- list(Ks, Linfs)
if (is.na(Linf_fix)) {
plot_dat <- reshape2::melt(score_mat)
#image(x = Linfs, y = Ks, z = t(score_mat), col = rsa.colors,
# ylab = "K", xlab = "Linf")
if (plot_title)
# title("Response surface analysis", line = 1)
if (contour) {
# contour(x = Linfs, y = Ks, z = t(score_mat), add = TRUE)
}
else {
if (is.numeric(contour)) {
#contour(x = Linfs, y = Ks, z = t(score_mat),
# add = TRUE, nlevels = contour)
}
else {
if (add.values) {
# text(x = plot_dat$Var2, y = plot_dat$Var1,
# round(as.numeric(plot_dat$value), digits = 2),
# cex = 0.6)
}
}
}
}
else {
if (all(Ks %in% exp(seq(log(0.1), log(10), length.out = 100)))) {
K_labels <- c(seq(0.1, 1, 0.1), 2:10)
K_plot <- log10(Ks)
K_ats <- log10(K_labels)
}
else {
K_labels <- Ks
K_plot <- Ks
K_ats <- Ks
}
phis <- round(log10(K_labels) + 2 * log10(Linfs), digits = 2)
op <- par(mar = c(12, 5, 4, 2))
# plot(K_plot, score_mat, type = "l", ylim = c(0, max(score_mat,
# na.rm = TRUE) * 1.4), ylab = "Score function", xlab = "Growth constant K (/year)",
# col = "red", lwd = 2, xaxt = "n")
#axis(1, at = K_ats, labels = K_labels)
#axis(1, at = K_ats, labels = phis, line = 5.5)
#mtext(text = expression(paste("Growth performance index (",
# phi, "')")), side = 1, line = 8.5)
#grid(nx = 0, NULL, lty = 6, col = "gray40")
#abline(v = K_ats, lty = 6, col = "gray40")
#if (plot_title)
# title("K-Scan", line = 1)
par(op)
}
Rn_max <- max(score_mat, na.rm = TRUE)[1]
idxs <- which(score_mat == Rn_max, arr.ind = TRUE)[1, ]
Linfest <- as.numeric(as.character(colnames(score_mat)[idxs[2]]))
Kest <- as.numeric(as.character(rownames(score_mat)[idxs[1]]))
tanchest <- as.numeric(as.character(ESP_tanch_L[idxs[1],
idxs[2]]))
phiL <- log10(Kest) + 2 * log10(Linfest)
pars <- list(Linf = Linfest, K = Kest, t_anchor = tanchest,
C = C, ts = ts, phiL = phiL)
final_res <- lfqFitCurves(lfq = res, par = pars, flagging.out = flagging.out,
agemax = agemax)
ret <- c(res, list(score_mat = score_mat, t_anchor_mat = ESP_tanch_L,
ncohort = final_res$ncohort, agemax = final_res$agemax,
par = pars, Rn_max = Rn_max))
class(ret) <- "lfq"
if (plot) {
#plot(ret, Fname = "rcounts")
Lt <- lfqFitCurves(ret, par = pars, draw = TRUE)
}
return(ret)
}
#
output <- RSA_plot_for_LINUX_no_plot_b(synLF5b, Linf_range = seq(78,88,0.2),
K_range = seq(0.3,0.7,0.1),
C = 0, ts = 0, MA = 9,
contour = 3, plot = FALSE )
output$par
output$Rn_max
outp.table[1,1] <- round(as.numeric(output$par[1]), 2)
outp.table[1,2] <- round(as.numeric(output$par[2]), 3)
outp.table[1,3] <- round(as.numeric(output$par[3]), 3)
outp.table[1,4] <- round(as.numeric(output$par[6]), 3)
outp.table[1,5] <- round(output$Rn_max, 4)
#setwd("~/WG_New_Length-Based_Methods")
write.table(outp.table,file = "exp5BOOT_1_outp_table_RSAcalbaMA9.csv" ,
append = TRUE,sep= ";", col.names = FALSE, row.names = FALSE)
# clear memory (avoids crashing the R session)
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
rm(list=ls()) # this will remove ALL objects (clear memory)
gc()
}
# setwd("~/WG_New_Length-Based_Methods")
boot.out51RSA <- read.table(file = "exp5BOOT_1_outp_table_RSAcalbaMA9.csv" , sep= ";")
names(boot.out51RSA) <-c("Linf",
"K",
"t_anchor",
"Phi_p" ,
"Rn_max" ,
"MA")
# Best fit parameters
boot.out11GA[which.max(boot.out11GA$Rn_max),]
# Conf.ints
quantile(boot.out11GA$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out11GA$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out11GA$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
################# RSA plot (1 plot, not for BOOTSTRAP)
# use one year only
synLF5b <- synLFQ5
synLF5b$catch <- synLF5b$catch[,1:12]
synLF5b$dates <- synLF5b$dates[1:12]
# For one-time use only
# Big picture and extremley slow RSA (wide range, high resolution) - not for BOOTSTRAP - just for one plot
#output.RSA.mat.for.plotSYN5b <- ELEFAN(synLF5b, Linf_range = seq(74,94,0.1),
# K_range = seq(0.3,0.7,0.1),
# C = 0, ts = 0,MA = 9,
# contour = 3, plot = F )
# "best fit" Rn values: par$Linf = 10.3, par$K = 1.4
saveRDS(output.RSA.mat.for.plotSYN5b,file = "output_RSA_mat_for_plotMA9syn5b.rds")
rm(output.RSA.mat.for.plot)
output.RSA.mat.for.plotSYN5b <- readRDS("output_RSA_mat_for_plotMA9syn5b.rds")
output.RSA.mat.for.plotSYN5b
#plot(output.RSA.mat.for.plotSYN5b)
# View(as.data.frame(output.RSA.mat.for.plotSYN5b$score_mat))
# plot(output.RSA.mat.for.plotSYN5b$score_mat)
##################
# FULL BOOTSTRAP 2 (narrow range, low resolution, no plot "fast" RSA plot)
#
# 10 runs take about 1 minute
for (i in 1:30)
{
# Sample and replace ALL monthly LFQ data
library(TropFishR)
# use one year only
synLF5b <- synLFQ5
synLF5b$catch <- synLF5b$catch[,1:12]
synLF5b$dates <- synLF5b$dates[1:12]
synLF5b.t <- synLF5b
ncol(synLF5b.t$catch) # n = x montlhy LFQ data
mlen <- synLF5b.t$midLengths # midlengths
c.int <- (mlen[2]-mlen[1]) #class interval (1 cm)
# "m" ist the month(i.e, the catch file column), 1 to 25
for (m in 1:(ncol(synLF5b.t$catch))) {
catch.m <- synLF5b.t$catch[,m]
#plot(catch.m ~ mlen, xlim = c(20,90),main = "original LFD")# original LFD
#rebuild original data
reb.catch.m <- rep(mlen,catch.m) # recreates a proxy of original data
# take a sample
sam.reb.catch.m<- sample(reb.catch.m,replace = TRUE)
hist.sam.reb.catch.m <- hist(sam.reb.catch.m, breaks = c(mlen-(c.int/2)), plot = FALSE)
#replaces data in LFQ file "2" with bootstrap sample data taken from "2"
synLF5b.t$catch[1:length(hist.sam.reb.catch.m$counts) ,m] <- hist.sam.reb.catch.m$counts
}
synLF5b.run1 <- synLF5b.t # new sample
RSA_plot_for_LINUX_no_plot_b <- function (x, Linf_fix = NA, Linf_range = NA, K_range = exp(seq(log(0.1),
log(10), length.out = 100)), C = 0, ts = 0,MA = 9, addl.sqrt = FALSE,
agemax = NULL, flagging.out = TRUE, method = "optimise",
cross.date = NULL, cross.midLength = NULL, cross.max = FALSE,
hide.progressbar = FALSE, plot = FALSE, contour = FALSE,
add.values = TRUE, rsa.colors = terrain.colors(20), plot_title = TRUE)
{
res <- x
classes <- res$midLengths
catch <- res$catch
dates <- res$dates
interval <- (classes[2] - classes[1])/2
n_samples <- dim(catch)[2]
n_classes <- length(classes)
if (is.na(Linf_fix) & is.na(Linf_range[1]))
Linf_range <- seq(classes[n_classes] - 5, classes[n_classes] +
5, 1)
res <- lfqRestructure(res, MA = MA, addl.sqrt = addl.sqrt)
catch_aAF_F <- res$rcounts
peaks_mat <- res$peaks_mat
ASP <- res$ASP
if (!is.na(Linf_fix)) {
Linfs <- Linf_fix
}
else Linfs <- Linf_range
Ks <- K_range
sofun <- function(tanch, lfq, par, agemax, flagging.out) {
Lt <- lfqFitCurves(lfq, par = list(Linf = par[1], K = par[2],
t_anchor = tanch, C = par[4], ts = par[5]), flagging.out = flagging.out,
agemax = agemax)
return(Lt$ESP)
}
ESP_tanch_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
ESP_list_L <- matrix(NA, nrow = length(Ks), ncol = length(Linfs))
writeLines(paste("Optimisation procuedure of ELEFAN is running. \nThis will take some time. \nThe process bar will inform you about the process of the calculations.",
sep = " "))
flush.console()
if (method == "cross") {
if (cross.max) {
max.rcount <- which.max(res$rcounts)
COMB <- expand.grid(midLengths = rev(rev(res$midLengths)),
dates = res$dates)
cross.date <- COMB$dates[max.rcount]
cross.midLength <- COMB$midLengths[max.rcount]
}
else {
if (is.null(cross.date) | is.null(cross.midLength)) {
stop("Must define both 'cross.date' and 'cross.midLength' when 'cross.max' equals FALSE")
}
}
}
if (!hide.progressbar) {
nlk <- prod(dim(ESP_list_L))
pb <- txtProgressBar(min = 1, max = nlk, style = 3)
counter <- 1
}
for (li in 1:length(Linfs)) {
ESP_tanch_K <- rep(NA, length(Ks))
ESP_list_K <- rep(NA, length(Ks))
for (ki in 1:length(Ks)) {
if (is.null(agemax)) {
agemax.i <- ceiling((1/-Ks[ki]) * log(1 - ((Linfs[li] *
0.95)/Linfs[li])))
}
else {
agemax.i <- agemax
}
if (method == "cross") {
t0s <- seq(floor(min(date2yeardec(res$dates)) -
agemax.i), ceiling(max(date2yeardec(res$dates))),
by = 0.01)
Ltlut <- Linfs[li] * (1 - exp(-(Ks[ki] * (date2yeardec(cross.date) -
t0s) + (((C * Ks[ki])/(2 * pi)) * sin(2 *
pi * (date2yeardec(cross.date) - ts))) - (((C *
Ks[ki])/(2 * pi)) * sin(2 * pi * (t0s - ts))))))
t_anchor <- t0s[which.min((Ltlut - cross.midLength)^2)]%%1
resis <- sofun(tanch = t_anchor, lfq = res,
par = c(Linfs[li], Ks[ki], NA, C, ts), agemax = agemax.i,
flagging.out = flagging.out)
ESP_list_K[ki] <- resis
ESP_tanch_K[ki] <- t_anchor
}
if (method == "optimise") {
resis <- optimise(f = sofun, lower = 0, upper = 1,
lfq = res, par = c(Linfs[li], Ks[ki], NA,
C, ts), agemax = agemax.i, flagging.out = flagging.out,
tol = 0.001, maximum = TRUE)
ESP_list_K[ki] <- resis$objective
ESP_tanch_K[ki] <- resis$maximum
}
if (!hide.progressbar) {
setTxtProgressBar(pb, counter)
counter <- counter + 1
}
}
ESP_list_L[, li] <- ESP_list_K
ESP_tanch_L[, li] <- ESP_tanch_K
}
dimnames(ESP_list_L) <- list(Ks, Linfs)
score_mat <- round((10^(ESP_list_L/ASP))/10, digits = 3)
rownames(score_mat) <- round(as.numeric(rownames(score_mat)),
digits = 2)
dimnames(ESP_tanch_L) <- list(Ks, Linfs)
if (is.na(Linf_fix)) {
plot_dat <- reshape2::melt(score_mat)
#image(x = Linfs, y = Ks, z = t(score_mat), col = rsa.colors,
# ylab = "K", xlab = "Linf")
if (plot_title)
# title("Response surface analysis", line = 1)
if (contour) {
# contour(x = Linfs, y = Ks, z = t(score_mat), add = TRUE)
}
else {
if (is.numeric(contour)) {
#contour(x = Linfs, y = Ks, z = t(score_mat),
# add = TRUE, nlevels = contour)
}
else {
if (add.values) {
# text(x = plot_dat$Var2, y = plot_dat$Var1,
# round(as.numeric(plot_dat$value), digits = 2),
# cex = 0.6)
}
}
}
}
else {
if (all(Ks %in% exp(seq(log(0.1), log(10), length.out = 100)))) {
K_labels <- c(seq(0.1, 1, 0.1), 2:10)
K_plot <- log10(Ks)
K_ats <- log10(K_labels)
}
else {
K_labels <- Ks
K_plot <- Ks
K_ats <- Ks
}
phis <- round(log10(K_labels) + 2 * log10(Linfs), digits = 2)
op <- par(mar = c(12, 5, 4, 2))
# plot(K_plot, score_mat, type = "l", ylim = c(0, max(score_mat,
# na.rm = TRUE) * 1.4), ylab = "Score function", xlab = "Growth constant K (/year)",
# col = "red", lwd = 2, xaxt = "n")
#axis(1, at = K_ats, labels = K_labels)
#axis(1, at = K_ats, labels = phis, line = 5.5)
#mtext(text = expression(paste("Growth performance index (",
# phi, "')")), side = 1, line = 8.5)
#grid(nx = 0, NULL, lty = 6, col = "gray40")
#abline(v = K_ats, lty = 6, col = "gray40")
#if (plot_title)
# title("K-Scan", line = 1)
par(op)
}
Rn_max <- max(score_mat, na.rm = TRUE)[1]
idxs <- which(score_mat == Rn_max, arr.ind = TRUE)[1, ]
Linfest <- as.numeric(as.character(colnames(score_mat)[idxs[2]]))
Kest <- as.numeric(as.character(rownames(score_mat)[idxs[1]]))
tanchest <- as.numeric(as.character(ESP_tanch_L[idxs[1],
idxs[2]]))
phiL <- log10(Kest) + 2 * log10(Linfest)
pars <- list(Linf = Linfest, K = Kest, t_anchor = tanchest,
C = C, ts = ts, phiL = phiL)
final_res <- lfqFitCurves(lfq = res, par = pars, flagging.out = flagging.out,
agemax = agemax)
ret <- c(res, list(score_mat = score_mat, t_anchor_mat = ESP_tanch_L,
ncohort = final_res$ncohort, agemax = final_res$agemax,
par = pars, Rn_max = Rn_max))
class(ret) <- "lfq"
if (plot) {
#plot(ret, Fname = "rcounts")
Lt <- lfqFitCurves(ret, par = pars, draw = TRUE)
}
return(ret)
}
# Big picture and veeery slow RSA (wide range, high resolution)
#output2.new <- RSA_plot_for_LINUX_no_plot_b(alba2.run1, Linf_range = seq(7,18,0.2),
# K_range = seq(0.4,7,0.2),
# C = 0, ts = 0, MA = 9,
# contour = 3, plot = FALSE )
# narrow and fast RSA (searches whithin narrow range)
output2.new <- RSA_plot_for_LINUX_no_plot_b(synLF5b.run1, Linf_range = seq(79,88,0.2),
K_range = seq(0.35,0.65,0.05),
C = 0, ts = 0, MA = 9,
contour = 3, plot = FALSE )
output2.new$par
output2.new$Rn_max
output2.new$par
#plot(output2.new)
#date.time <- date()
Linf_out <- round(as.numeric(output2.new$par[1]), 2)
K_out <- round(as.numeric(output2.new$par[2]), 3)
t_anchor <- round(as.numeric(output2.new$par[3]), 3)
Phi_p <- round(as.numeric(output2.new$par[6]), 3)
Rn_max <- round(as.numeric(output2.new$Rn_max),3)
output.table <- data.frame(Linf_out, K_out, t_anchor, Phi_p,Rn_max)
write.table(output.table, file = "exp6BOOT_2__outp_table_RSAsynlf5b_MA9.csv",
sep = ";",col.names = FALSE, row.names = FALSE, append = TRUE )
rm(list=ls())
gc()
} # end of loop
boot.out62RSA_syn5b <- read.table( file = "exp6BOOT_2__outp_table_RSAsynlf5b_MA9.csv", sep = ";")
names(boot.out62RSA_syn5b) <- c("Linf", "K", "t_anchor", "Phi_p","Rn_max")
boot.out62RSA_syn5b
# Best fit parameters
boot.out62RSA_syn5b[which.max(boot.out62RSA_syn5b$Rn_max),]
# Conf.ints
quantile(boot.out62RSA_syn5b$K ,c(0.025, 0.975)) # OK , _SA: 95 conf. int K = to
quantile(boot.out62RSA_syn5b$Linf ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Linf = to
quantile(boot.out62RSA_syn5b$Phi_p ,c(0.025, 0.975)) # OK , _SA: 95 conf. int Phi_Prime = to
##
### END OF SCRIPT ###
Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.