R/plot_eco_hr_act_eco_optimization_diagnostics.R

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

#
# plot_eco_hr_act_eco_optimization_diagnostics.R
#
#' Create diagnostic plots from parameter optimization procedures. 
#'
#' Create a plot showing the distributions of parameter proposals relative to the initial values in an optimization procedure, and the final accepted value relative to the initial. 
#' The function accepts data from the optimization function versions for optimizing ecology model parameters, harvest ratio parameters, and activity rate parameters (the latter against ecosystem state targets).
#'
#' In the case of the plot derived from ecology parameter optimization, the parameters are identified by an abbreviated name. A key to the paramaters identities
#' can be downloaded in dataframe format by the function e2e_get_parmdoc().
#' 
#' @param model R-list object defining the model configuration used to generate the data and compiled by the e2e_read() function.
#' @param selection Text string from a list identifying which type of optimization procedure generated the data to be plotted. Select from: "ECO", "HR", "ACT", corresponding to the functions e2e_optimize_eco(), e2e_optimize_hr(), and e2e_optimize_act(). Remember to include the phrase within "" quotes.
#' @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 Graphical display in a new graphics window, and a list object of the plotted data.
#'
#' @noRd
#
# ------------------------------------------------------------------------------

plot_eco_hr_act_eco_optimization_diagnostics <- function(model, selection="ECO",use.saved=FALSE, use.example=FALSE, results=NULL) {

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

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

	gear_labels     <- elt(model, "data", "fleet.model", "gear_labels")
	gear_codes     <- elt(model, "data", "fleet.model", "gear_codes")


allowed.list<-c("ECO","HR","ACT")
if(selection %in% allowed.list==FALSE){
	message("The argument 'selection' must be one of the following...")
	message("  'ECO' = results from optimizing ecology parameters")
	message("  'HR'  = results from optimizing harvest ratio multipliers")
	message("  'ACT' = results from optimizing gear activity multipliers (against observed ecosystem data)")
	stop(" unknown selection '", selection, "' !\n")
}

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 & selection=="ACT" ){
        stop("No example data available for activity fitting to ecosystem targets ! \n")
}

if(use.example==TRUE & selection=="ECO" & (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 & selection=="HR" & (model.variant=="2003-2013")==FALSE ){
        stop("Example data only available for the 2003-2013 variant of the North Sea model ! \n")
}

if(selection=="ECO"){
	coreproposalname<-"annealing_par_proposalhistory"
	coreacceptedname<-"annealing_par_acceptedhistory"
}

if(selection=="HR"){
	coreproposalname<-"annealing_HRmult_proposalhistory"
	coreacceptedname<-"annealing_HRmult_acceptedhistory"
	parmlist<-c("Planktivorous fish","Demersal fish",
		     "Migratory fish","Susp.dep. benthos",
		     "Carn/scav. benthos","Carniv.zooplankton",
		     "Birds","Pinnipeds","Cetaceans","Macrophytes")
}

if(selection=="ACT"){
	coreproposalname<-"annealing_ACTmult_proposalhistory"
	coreacceptedname<-"annealing_ACTmult_acceptedhistory"
	parmlist<-gear_labels
#	parmlist<-gear_codes
}

if(use.example==TRUE){
	proposaldata <- get.example.results(model.name, model.variant, coreproposalname)
	accepteddata <- get.example.results(model.name, model.variant, coreacceptedname)
}

if(use.saved==TRUE){
	proposalfile	<- csvname(resultsdir, coreproposalname, model.ident)
	acceptedfile	<- csvname(resultsdir, coreacceptedname, model.ident)
	message("Reading saved data from folder ", resultsdir)

	proposaldata<- readcsv(proposalfile)
	accepteddata<- readcsv(acceptedfile)
}

if(is.list(results)==TRUE & use.saved==FALSE & use.example==FALSE){
	proposaldata <- results$parameter_proposal_history
	accepteddata <- results$parameter_accepted_history
}


	Nind<-ncol(proposaldata)-1
	Niter<-nrow(proposaldata)


#extract the initial and final accepted parameter values

initial.values <-as.numeric(accepteddata[1,1:Nind])
final.values   <-as.numeric(accepteddata[Niter,1:Nind])

#Calculate proposals relative to initial values
proposaldata_dif<-proposaldata[,1:Nind]
for(j in 1:Nind){
proposaldata_dif[,j] <- (proposaldata[,j]-initial.values[j])/initial.values[j]
}

final_dif<-(final.values-initial.values)/initial.values

	if(selection=="ECO") parmlist<-colnames(proposaldata_dif)

	FINAL_array<-array(dim=c(5,Nind),NA)
	colnames(FINAL_array)<-parmlist
	PROPOSAL_array<-FINAL_array
	rownames(PROPOSAL_array)<-c("0.005","0.25","median","0.75","0.995")
	rownames(FINAL_array)<-c("maxlikvalue","maxlikvalue","maxlikvalue","maxlikvalue","maxlikvalue")

	for(j in 1:Nind){
	if(is.na(final_dif[j])==FALSE){
	FINAL_array[,j]<-final_dif[j]
	PROPOSAL_array[,j]<-quantile(proposaldata_dif[,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_PROP<-list(stats=PROPOSAL_array,conf=NULL,out=numeric(length=0),names=parmlist)
	bxpdata_FINAL<-list(stats=FINAL_array,conf=NULL,out=numeric(length=0),names=parmlist)

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

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

	par(mfrow=c(1,1))
	if(selection=="ECO") par(mar=c(3.2,8,1,1))
	if(selection=="HR") par(mar=c(3.2,10,1,1))
	if(selection=="ACT") par(mar=c(3.2,14,1,1))
  	bxp(bxpdata_PROP,boxwex=0.75,at=seq(0,Nind-1),horizontal=TRUE,show.names=FALSE,las=1,ylim=c(-axmax,axmax),boxcol="black",whiskcol="black",whisklty="solid",medcol="black",staplecol="black")
  	if(selection=="ECO") bxp(bxpdata_FINAL,boxwex=1.25,add=TRUE,at=seq(0,Nind-1),horizontal=TRUE,show.names=FALSE,las=1,boxcol="red",whiskcol="red",whisklty="solid",medcol="red",staplecol="red")
  	if(selection=="HR" | selection=="ACT") bxp(bxpdata_FINAL,boxwex=1,add=TRUE,at=seq(0,Nind-1),horizontal=TRUE,show.names=FALSE,las=1,boxcol="red",whiskcol="red",whisklty="solid",medcol="red",staplecol="red")
	if(selection=="ECO") axis(parmlist,side=2,las=1,at=seq(0,Nind-1),cex.axis=0.38)
	if(selection=="HR") axis(parmlist,side=2,las=1,at=seq(0,Nind-1),cex.axis=1)
	if(selection=="ACT") axis(parmlist,side=2,las=1,at=seq(0,Nind-1),cex.axis=0.8)
	mtext("Parameters relative to initial values",cex=1.1,side=1,line=2)
	abline(v=0,lty="dashed")

	out_data<-list(proposed_value_disribution_rel_initial=PROPOSAL_array,finalvalue_rel_initial=FINAL_array[1,])

return(out_data)

}