Nothing
#
# 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.