R/plot_act_hr_optimization_diagnostics.R

In StrathE2E2: End-to-End Marine Food Web Model

#
# plot_act_hr_optimization_diagnostics.R
#
#' Create a diagnostic plot to show the performance of the simulated annealing function to optimize gear activity rates.
#'
#' Simulated annealing is an iterative random-walk type process which searches the parameter space
#' of a model to locate the combination which maximises the likelihood of a set of observed
#' data corresponding to a suite of derived outputs. Parameter combinations which result in an improved likelihood
#' may be rejected according to a probability ('temperature') which decreases as the iterations progress. This is to  
#' avoid becoming stuck at local likelihood-maxima. The rate at which the 'temperature' decreases is set
#' by a 'cooling' parameter (fraction of previous temperature at each iteration, 0<value<1). The task here is a bit harder than normal
#' because of the potentially large overlap in the selectivity patterns of the fishing gears with respect to the living guilds in the ecosystem. So there are
#' some non-standard features to try and avoid local maxima.
#'
#' The function produces a diagnostic box and whisker plots showing the uncertainty in the estimated activity rate multipliers
#' and in the estimated harvest ratios across the end-points of all the trajectories. In each diagnostic plot the maximum likelihood model resullts are shown by a red
#' bar, and the distribution of all results by black box and whiskers.
#'
#' @param model R-list object defining the model configuration used to generate the data and compiled by the e2e_read() function.
#' @param use.saved Logical. If TRUE use data from a prior user-defined run held as csv files data in the current results folder, (default=FALSE).
#' @param use.example Logical. If TRUE use pre-computed example data from the internal North Sea model rather than user-generated data, (default=FALSE).
#' @param results List object generated by the function e2e_optimize_act() with the argument selection="HR", (default = NULL)
#'
#' @return List object containing processed data and graphical display of diagnostic data
#'
#' @importFrom graphics axis legend lines mtext par plot points boxplot
#' @importFrom stats quantile
#'
#' @noRd
#
# ------------------------------------------------------------------------------

plot_act_hr_optimization_diagnostics <- function(model, use.saved=FALSE, use.example=FALSE, results=NULL) {

start_par = par()$mfrow
on.exit(par(mfrow = start_par))

	model.name 	<- elt(model, "setup", "model.name")
	model.variant 	<- elt(model, "setup", "model.variant")


if(use.saved==TRUE & use.example==TRUE){
        stop("use.saved and use.example cannot both be TRUE ! \n")
}

if(use.saved==TRUE & is.list(results)==TRUE){
        stop("use.saved is TRUE but a dataframe object has also been specified ! \n")
}

if(use.example==TRUE & is.list(results)==TRUE){
        stop("use.example is TRUE but a dataframe object has also been specified ! \n")
}

if(use.saved==FALSE & use.example==FALSE & is.list(results)==FALSE){
        stop("no source of data has been specified ! \n")
}

if(use.example==TRUE & (model.name=="North_Sea")==FALSE ){
        stop("Example data only available for the North Sea model ! \n")
}

if(use.example==TRUE & (model.variant=="1970-1999")==FALSE ){
        stop("Example data only available for the 1970-1999 variant of the North Sea model ! \n")
}


if(use.example==TRUE){
	model.name 	<- elt(model, "setup", "model.name")
	model.variant 	<- elt(model, "setup", "model.variant")
	gearmultrelative_results <- get.example.results(model.name, model.variant, "activity_optim_gearmult_relinitial_history")
	HRvaluesrelative_results <- get.example.results(model.name, model.variant, "activity_optim_harvestratio_reltarget_history")
}

if(use.saved==TRUE) {
	resultsdir	<- elt(model, "setup", "resultsdir")
	model.ident	<- elt(model, "setup", "model.ident")
	hrdatafile   <- csvname(resultsdir, "activity_optim_harvestratio_reltarget_history", model.ident)
	multdatafile <- csvname(resultsdir, "activity_optim_gearmult_relinitial_history", model.ident)
	check.exists(hrdatafile)
	check.exists(multdatafile)
	print(paste("Using data held in files with model.ident ",model.ident," from a past model run"))
	gearmultrelative_results<-readcsv(multdatafile)
	HRvaluesrelative_results<-readcsv(hrdatafile)
}

if(use.saved==FALSE & use.example==FALSE & is.list(results)==TRUE){
	gearmultrelative_results   <- results$gear_mult_rel_initial
	HRvaluesrelative_results <- results$harvest_ratio_rel_target
}

#--------------------------------------------------------

	# Now, find the best fit values of gearmultrelative and HRrelative
	gearmultrelative_BEST<-gearmultrelative_results[which(gearmultrelative_results[,13]==max(gearmultrelative_results[,13])),]
	HRvaluesrelative_BEST<-HRvaluesrelative_results[which(HRvaluesrelative_results[,21]==max(HRvaluesrelative_results[,21])),]

	#Set up an array with 5 rows to hold the BEST data - the data in each row is going to be the same
	# - this is just a trick to make plotting a bit easier later
	gearmult_BEST_array_relative<-array(dim=c(5,12),NA)
	HRvalues_BEST_array_relative<-array(dim=c(5,20),NA)
	gear_code_list<-colnames(gearmultrelative_results)
	guild_list<-colnames(HRvaluesrelative_results)
	colnames(gearmult_BEST_array_relative)<-gear_code_list[1:12] 
	colnames(HRvalues_BEST_array_relative)<-guild_list[1:20]
	rownames(gearmult_BEST_array_relative)<-c("maxlikvalue","maxlikvalue","maxlikvalue","maxlikvalue","maxlikvalue")
	rownames(HRvalues_BEST_array_relative)<-c("maxlikvalue","maxlikvalue","maxlikvalue","maxlikvalue","maxlikvalue")

	for(j in 1:12){
		gearmult_BEST_array_relative[,j]<-rep(as.numeric(gearmultrelative_BEST[j]),5)
	}

	for(j in 1:20){
		HRvalues_BEST_array_relative[,j]<-rep(as.numeric(HRvaluesrelative_BEST[j]),5)
	}


	#Set up an array with 5 rows to hold the quantiles of the history data
	gearmult_CREDINT_array_relative<-array(dim=c(5,12),NA)
	HRvalues_CREDINT_array_relative<-array(dim=c(5,20),NA)
	colnames(gearmult_CREDINT_array_relative)<-gear_code_list[1:12] 
	colnames(HRvalues_CREDINT_array_relative)<-guild_list[1:20]
	rownames(gearmult_CREDINT_array_relative)<-c("0.005","0.25","median","0.75","0.995")
	rownames(HRvalues_CREDINT_array_relative)<-c("0.005","0.25","median","0.75","0.995")

	for(j in 1:12){
	if(is.na(gearmultrelative_BEST[j])==FALSE){
	gearmult_CREDINT_array_relative[,j]<-quantile(gearmultrelative_results[,j],probs=c(0.005,0.25,0.5,0.75,0.995))
	}
	}

	for(j in 1:20){
	if(is.na(HRvaluesrelative_BEST[j])==FALSE){
	HRvalues_CREDINT_array_relative[,j]<-quantile(HRvaluesrelative_results[,j],probs=c(0.005,0.25,0.5,0.75,0.995))
	}
	}

	#Set up the list objects needed by the bxp function which draws box and whisker plots

	bxpdata_G<-list(stats=gearmult_CREDINT_array_relative,conf=NULL,out=numeric(length=0),names=colnames(gearmult_CREDINT_array_relative))
	bxpdata_Gbest<-list(stats=gearmult_BEST_array_relative,conf=NULL,out=numeric(length=0),names=colnames(gearmult_CREDINT_array_relative))


	bxpdata_H<-list(stats=HRvalues_CREDINT_array_relative,conf=NULL,out=numeric(length=0),names=colnames(HRvalues_CREDINT_array_relative))
	bxpdata_Hbest<-list(stats=HRvalues_BEST_array_relative,conf=NULL,out=numeric(length=0),names=colnames(HRvalues_CREDINT_array_relative))

	#-------------------------------------------------------

	#Plot the gear multipler results 


	par(mfrow=c(1,2))

	plymax<-max(gearmult_CREDINT_array_relative,na.rm=TRUE)
	plymin<-min(gearmult_CREDINT_array_relative,na.rm=TRUE)

	axmax<-sqrt(max(c(plymax^2,plymin^2)))

	par(mar=c(3.2,3,1,2.5))
  	bxp(bxpdata_G,boxwex=0.35,at=seq(0,11),horizontal=TRUE,show.names=FALSE,las=1,ylim=c(-axmax,axmax),boxcol="black",whiskcol="black",whisklty="solid",medcol="black",staplecol="black")
  	bxp(bxpdata_Gbest,boxwex=0.35,add=TRUE,at=seq(0,11),horizontal=TRUE,show.names=FALSE,las=1,boxcol="red",whiskcol="red",whisklty="solid",medcol="red",staplecol="red")
	axis(colnames(gearmult_CREDINT_array_relative),side=2,las=1,at=seq(0,11),cex.axis=1)
	mtext("Gear act. multipliers relative to initial",cex=1.1,side=1,line=2)
	abline(v=0,lty="dashed")

	#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


	#Plot the Harvest ratio results 

	plymax<-max(HRvalues_CREDINT_array_relative,na.rm=TRUE)
	plymin<-min(HRvalues_CREDINT_array_relative,na.rm=TRUE)

	axmax<-sqrt(max(c(plymax^2,plymin^2)))

	par(mar=c(3.2,4,1,1))
  	bxp(bxpdata_H,boxwex=0.35,at=seq(0,19),horizontal=TRUE,show.names=FALSE,las=1,ylim=c(-axmax,axmax),boxcol="black",whiskcol="black",whisklty="solid",medcol="black",staplecol="black")
  	bxp(bxpdata_Hbest,boxwex=0.35,add=TRUE,at=seq(0,19),horizontal=TRUE,show.names=FALSE,las=1,boxcol="red",whiskcol="red",whisklty="solid",medcol="red",staplecol="red")
	axis(colnames(HRvalues_CREDINT_array_relative),side=2,las=1,at=seq(0,19),cex.axis=0.8)
	mtext("Harvest ratios relative to target",cex=1.1,side=1,line=2)
	abline(v=0,lty="dashed")
	

	#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

output_data <- list(gearmult_distrubution_rel_initial = gearmult_CREDINT_array_relative,
		    gearmult_maxlikelihood_rel_initial= gearmult_BEST_array_relative[1,],
		    HRvalues_distrubution_rel_initial = HRvalues_CREDINT_array_relative,
		    HRvalues_maxlikelihood_rel_initial= HRvalues_BEST_array_relative[1,])

return(output_data)

}