#' @title ELEFAN_SA
#'
#' @description Electronic LEngth Frequency ANalysis with simulated annealing
#' for estimating growth parameters.
#'
#' @param lfq a list consisting of following parameters:
#' \itemize{
#' \item \strong{midLengths} midpoints of the length classes,
#' \item \strong{dates} dates of sampling times (class Date),
#' \item \strong{catch} matrix with catches/counts per length class (row)
#' and sampling date (column);
#' }
#' @param seasonalised logical; indicating if the seasonalised von Bertalanffy
#' growth function should be applied (default: FALSE).
#' @param init_par a list providing the Initial values for the components to be
#' optimized. When set to NULL the following default values are used:
#' \itemize{
#' \item \strong{Linf} length infinity in cm (default is the maximum
#' length class in the data),
#' \item \strong{K} curving coefficient (default: 0.5),
#' \item \strong{t_anchor} time point anchoring growth curves in year-length
#' coordinate system, corrsponds to peak spawning month (range: 0 to 1, default: 0.5),
#' \item \strong{C} amplitude of growth oscillation (range: 0 to 1, default: 0),
#' \item \strong{ts} summer point (ts = WP - 0.5) (range: 0 to 1, default: 0);
#' }
#' @param low_par a list providing the lower bounds for components. When set to
#' NULL the following default values are used:
#' \itemize{
#' \item \strong{Linf} length infinity in cm (default is calculated from maximum
#' length class in the data),
#' \item \strong{K} curving coefficient (default: 0.01),
#' \item \strong{t_anchor} time point anchoring growth curves in year-length
#' coordinate system, corrsponds to peak spawning month (range: 0 to 1, default: 0),
#' \item \strong{C} amplitude of growth oscillation (range: 0 to 1, default: 0),
#' \item \strong{ts} summer point (ts = WP - 0.5) (range: 0 to 1, default: 0);
#' }
#' @param up_par a list providing the upper bounds for components. When set to
#' NULL the following default values are used:
#' \itemize{
#' \item \strong{Linf} length infinity in cm (default is calculated from maximum
#' length class in the data),
#' \item \strong{K} curving coefficient (default: 0.01),
#' \item \strong{t_anchor} time point anchoring growth curves in year-length
#' coordinate system, corrsponds to peak spawning month (range: 0 to 1, default: 0),
#' \item \strong{C} amplitude of growth oscillation (range: 0 to 1, default: 0),
#' \item \strong{ts} summer point (ts = WP - 0.5) (range: 0 to 1, default: 0);
#' }
#' @param SA_time numeric; Maximum running time in seconds (default : 60 * 1).
#' @param maxit Integer. Maximum number of iterations of the
#' algorithm. Default is NULL.
#' @param nb.stop.improvement Integer. The program will stop when
#' there is no any improvement in 'nb.stop.improvement'
#' steps. Default is NULL
#' @param SA_temp numeric; Initial value for temperature (default : 1e5).
#' @param verbose logical; TRUE means that messages from the algorithm
#' are shown (default : TRUE).
#' @param MA number indicating over how many length classes the moving average
#' should be performed (defalut: 5, for
#' more information see \link{lfqRestructure}).
#' @param addl.sqrt Passed to \link{lfqRestructure}. Applied an additional
#' square-root transformation of positive values according to Brey et al. (1988).
#' (default: FALSE, for more information see \link{lfqRestructure}).
#' @param agemax maximum age of species; default NULL, then estimated from Linf
#' @param flagging.out logical; passed to \link{lfqFitCurves}. Default is TRUE
#' @param plot logical; Plot restructured counts with fitted lines using
#' \code{\link{plot.lfq}} and \code{\link{lfqFitCurves}} (default : FALSE).
#' @param plot.score logical; Plot simulated annealing score progression.
#' (Default: plot.score=TRUE)
#'
#' @examples
#' \donttest{
#' ## synthetic lfq data example
#' data(synLFQ4)
#' plot(synLFQ4, Fname="catch")
#'
#' # ELEFAN_SA (takes approximately 2 minutes)
#' output <- ELEFAN_SA(synLFQ4, SA_time = 60*2, seasonalised = TRUE, MA = 11,
#' init_par = list(Linf = 75, K = 0.5, t_anchor = 0.5, C = 0.5, ts = 0.5),
#' low_par = list(Linf = 70, K = 0.3, t_anchor = 0, C = 0, ts = 0),
#' up_par = list(Linf = 90, K = 0.7, t_anchor = 1, C = 1, ts = 1))
#' output$par
#' output$Rn_max
#'
#' # view fit
#' plot(output)
#'
#' # or
#' plot(output, draw = FALSE)
#' lfqFitCurves(output, col=1, par=output$par, draw=TRUE)$ESP
#'
#' # compare to original parameters
#' tmp <- lfqFitCurves(output, col=4, lty=1,
#' par=list(Linf=80, K=0.5, t_anchor=0.25, C=0.75, ts=0.5), draw=TRUE)
#' tmp$fESP
#' output$Rn_max
#' }
#'
#' @details A more detailed description of the simulated annealing (SA) can be found in
#' Xiang et al. (2013). The score value \code{cost_value} is not comparable with
#' the score value of the other ELEFAN functions (\code{\link{ELEFAN}} or
#' \code{\link{ELEFAN_GA}}).
#'
#' @return A list with the input parameters and following list objects:
#' \itemize{
#' \item \strong{rcounts}: restructured frequencies,
#' \item \strong{peaks_mat}: matrix with positive peaks with distinct values,
#' \item \strong{ASP}: available sum of peaks, sum of posititve peaks which
#' could be potential be hit by
#' growth curves,
#' \item \strong{ncohort}: maximum age of species,
#' \item \strong{agemax}: maximum age of species,
#' \item \strong{par}: a list with the parameters of the von Bertalanffy growth
#' function:
#' \itemize{
#' \item \strong{Linf}: length infinity in cm,
#' \item \strong{K}: curving coefficient;
#' \item \strong{t_anchor}: time point anchoring growth curves in year-length
#' coordinate system, corrsponds to peak spawning month,
#' \item \strong{C}: amplitude of growth oscillation
#' (if \code{seasonalised} = TRUE),
#' \item \strong{ts}: summer point of oscillation (ts = WP - 0.5)
#' (if \code{seasonalised} = TRUE),
#' \item \strong{phiL}: growth performance index defined as
#' phiL = log10(K) + 2 * log10(Linf);
#' }
#' \item \strong{Rn_max}: highest score value (absolute value of cost function,
#' comparable with ELEFAN and ELEFAN_GA).
#' }
#'
#' @importFrom graphics par plot title lines grid
#' @importFrom grDevices adjustcolor
#' @importFrom stats median
#' @importFrom GenSA GenSA
#' @importFrom utils flush.console
#'
#' @references
#' Brey, T., Soriano, M., and Pauly, D. 1988. Electronic length frequency analysis: a revised and
#' expanded
#' user's guide to ELEFAN 0, 1 and 2.
#'
#' Pauly, D. and N. David, 1981. ELEFAN I, a BASIC program for the objective extraction of
#' growth parameters from length-frequency data. \emph{Meeresforschung}, 28(4):205-211
#'
#' Xiang, Y., Gubian, S., Suomela, B., & Hoeng, J. (2013). Generalized simulated
#' annealing for global optimization: the GenSA Package. R Journal, 5(1), 13-28.
#' @export
ELEFAN_SA <- function(lfq,
seasonalised = FALSE,
init_par = list(Linf = 50, K = 0.5, t_anchor = 0.5, C = 0, ts = 0),
low_par = NULL,
up_par = NULL,
SA_time = 60 * 1,
maxit = NULL,
nb.stop.improvement = NULL,
SA_temp = 1e5,
verbose = TRUE,
MA = 5, addl.sqrt = FALSE,
agemax = NULL,
flagging.out = TRUE,
plot = FALSE,
plot.score = TRUE){
res <- lfq
classes <- res$midLengths
n_classes <- length(classes)
Linf_est <- classes[n_classes]
# initial parameters
init_par_ALL <- list(Linf = Linf_est,
K = 0.5, t_anchor = 0.5, C = 0, ts = 0)
init_Linf <- ifelse("Linf" %in% names(init_par),
get("Linf", init_par),
get("Linf", init_par_ALL))
init_K <- ifelse("K" %in% names(init_par),
get("K", init_par),
get("K", init_par_ALL))
init_tanc <- ifelse("t_anchor" %in% names(init_par),
get("t_anchor", init_par),
get("t_anchor", init_par_ALL))
init_C <- ifelse("C" %in% names(init_par),
get("C", init_par),
get("C", init_par_ALL))
init_ts <- ifelse("ts" %in% names(init_par),
get("ts", init_par),
get("ts", init_par_ALL))
# lower parameter bounds
low_par_ALL <- list(Linf = Linf_est * 0.5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0)
low_Linf <- ifelse("Linf" %in% names(low_par),
get("Linf", low_par),
get("Linf", low_par_ALL))
low_K <- ifelse("K" %in% names(low_par),
get("K", low_par),
get("K", low_par_ALL))
low_tanc <- ifelse("t_anchor" %in% names(low_par),
get("t_anchor", low_par),
get("t_anchor", low_par_ALL))
low_C <- ifelse("C" %in% names(low_par),
get("C", low_par),
get("C", low_par_ALL))
low_ts <- ifelse("ts" %in% names(low_par),
get("ts", low_par),
get("ts", low_par_ALL))
# upper parameter bounds
up_par_ALL <- list(Linf = Linf_est * 1.5,
K = 1,
t_anchor = 1,
C = 1,
ts = 1)
up_Linf <- ifelse("Linf" %in% names(up_par),
get("Linf", up_par),
get("Linf", up_par_ALL))
up_K <- ifelse("K" %in% names(up_par),
get("K", up_par),
get("K", up_par_ALL))
up_tanc <- ifelse("t_anchor" %in% names(up_par),
get("t_anchor", up_par),
get("t_anchor", up_par_ALL))
up_C <- ifelse("C" %in% names(up_par),
get("C", up_par),
get("C", up_par_ALL))
up_ts <- ifelse("ts" %in% names(up_par),
get("ts", up_par),
get("ts", up_par_ALL))
# ELEFAN 0
res <- lfqRestructure(res, MA = MA, addl.sqrt = addl.sqrt)
catch_aAF_F <- res$rcounts
peaks_mat <- res$peaks_mat
ASP <- res$ASP
# seasonalised cost function
soSAfun <- function(lfq, par=c(init_Linf, init_K, init_tanc, init_C, init_ts),
agemax, flagging.out){
Lt <- lfqFitCurves(lfq,
par=list(Linf=par[1], K=par[2], t_anchor=par[3], C=par[4], ts=par[5]),
flagging.out = flagging.out, agemax = agemax)
return(-Lt$fESP)
}
# cost function
SAfun <- function(lfq, par=c(init_Linf, init_K, init_tanc),
agemax, flagging.out){
Lt <- lfqFitCurves(lfq,
par=list(Linf=par[1], K=par[2], t_anchor=par[3], C = 0, ts = 0),
flagging.out = flagging.out, agemax = agemax)
return(-Lt$fESP)
}
## control list
control <- list(temperature = SA_temp,
verbose = verbose)
if(!is.null(SA_time)) control$max.time = SA_time
if(!is.null(maxit)) control$maxit = maxit
if(!is.null(nb.stop.improvement)) control$nb.stop.improvement
if(seasonalised){
# Simulated annealing with seasonalised VBGF
writeLines(paste(
"Simulated annealing is running. \nThis will take approximately",
round(SA_time/60,digits=2),
"minutes.",
sep=" "
))
flush.console()
SAfit <- GenSA::GenSA(
par = c(init_Linf, init_K, init_tanc, init_C, init_ts),
fn = soSAfun,
lower = c(low_Linf, low_K, low_tanc, low_C, low_ts),
upper = c(up_Linf, up_K, up_tanc, up_C, up_ts),
agemax = agemax,
flagging.out = flagging.out,
control = control,
lfq = res
)
pars <- as.list(SAfit$par)
names(pars) <- c("Linf","K","t_anchor", "C", "ts")
}else{
# Simulated annealing
writeLines(paste(
"Simulated annealing is running. \nThis will take approximately",
round(SA_time/60,digits=2),
"minutes.",
sep=" "
))
flush.console()
SAfit <- GenSA::GenSA(par = c(init_Linf, init_K, init_tanc),
fn = SAfun,
lower = c(low_Linf, low_K, low_tanc),
upper = c(up_Linf, up_K, up_tanc),
agemax = agemax,
flagging.out = flagging.out,
control = control,
lfq = res
)
pars <- as.list(SAfit$par)
names(pars) <- c("Linf","K","t_anchor")
}
# Score graph in GA style
## CHECK: trace.mat not reported any longer from GenSA, written to file but nr == 0 ...
if(FALSE){
tmp <- as.data.frame(SAfit$trace.mat)
meani <- aggregate(tmp$function.value, list(step = tmp$nb.steps),mean, na.rm = TRUE)
exe <- aggregate(tmp$current.minimum, list(step = tmp$nb.steps),mean, na.rm = TRUE)
medi <- aggregate(tmp$function.value, list(step = tmp$nb.steps),median, na.rm = TRUE)
ylim <- c(min(range(exe$x,na.rm = TRUE, finite = TRUE)),
max(range(meani$x, na.rm = TRUE, finite = TRUE)))
if(plot.score){
op <- par(mar=c(5.1, 4.1, 1, 4.1))
plot(tmp$nb.steps, tmp$function.value, type = "n", ylim = ylim, xlab = "Iteration",
ylab = "Cost value")
graphics::grid(equilogs = FALSE)
points(tmp$nb.steps, tmp$current.minimum, type = "o", pch = 16, lty = 1,
col = "green3", cex = 0.7)
points(meani$step, meani$x, type = "o", pch = 1, lty = 2,
col = "dodgerblue3", cex = 0.7)
polygon(c(meani$step, rev(meani$step)),
c(exe$x, rev(medi$x)),
border = FALSE, col = adjustcolor("green3", alpha.f = 0.1))
par(new=TRUE)
plot(tmp$nb.steps, tmp$temperature, t="l", col=2, lty=2, log="y", axes = FALSE, xlab = "", ylab = "")
axis(4, col=2, col.axis=2); mtext(text = "Temperature", side = 4, line = par()$mgp[1], col=2)
legend("topright", legend = c("Best", "Mean", "Median", "Temperature"),
col = c("green3", "dodgerblue3", adjustcolor("green3", alpha.f = 0.1), 2),
pch = c(16, 1, NA, NA), lty = c(1,2,1,2),
lwd = c(1, 1, 10, 1), pt.cex = c(rep(0.7,2), 2, NA),
inset = 0.02)
par(op)
}
}
final_res <- lfqFitCurves(lfq = res,par=pars,flagging.out = flagging.out,
agemax = agemax)
# growth performance index
phiL <- log10(pars$K) + 2 * log10(pars$Linf)
pars$phiL <- phiL
# Results
res$ncohort <- final_res$ncohort
res$agemax <- final_res$agemax
res$par <- pars
res$fESP <- abs(SAfit$value)
res$Rn_max <- abs(SAfit$value)
if(plot){
plot(res, Fname = "rcounts")
Lt <- lfqFitCurves(res, par = res$par, draw=TRUE)
}
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.