R/e2e_run_sens.R

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

#
# e2e_run_sens.R
#'
#' Perform a parameter sensitivity analysis on the StrathE2E model.
#'
#' Performs a one-at-a-time parameter sensitivity analysis on the StrathE2E model using
#' the Morris Method for factorial sampling of the physical configuration parameters, ecology model
#' parameters, the fishing fleet model parameters, and the
#' environmental forcings.
#'
#' The basis for the method is a scheme for sampling the model parameters, one at a time, from distributions around baseline sets, and testing the effects on the performance
#' of the model against some criterion. The default criterion is the likelihood of the observational data on the state of the ecosystem that
#' are used as the target for parameter optimization in the various simulated annealing functions supplied with the package. However, the criterion can in principle 
#' be any metric output from the model, e.g. a state variable value at some point in time or averaged over a time interval, or one of the computed fluxes.
#'
#' The process requires an initial set of parameters for the model. We refer to this as the 'parent' parameter set. It is recommended that this should be the parameters producing
#' the maximum likelihood of the observational target data (as estimated by e.g. the e2e_optimize_eco() function). The MODEL_SETUP.csv file in the folder
#' /Models/Modelname/Modelvariant/ should be configured to point to the relevant files, and then loaded with the e2e_read() function. 
#'
#' From this parent set, a series of 'child' parameter sets are generated by applying a separate random increment to each parameter drawn from a
#' gaussian distribution of mean 0 and standard deviation given by a fixed coefficient
#' of variation applied to the parent-set value of each parameter. 
#'
#' For each of the child-sets of parameters, the chosen output criterion is saved following runs of StrathE2E to stationary state.
#' We refer to these as trajectory baselines.
#'
#' Then, for each trajectory, the parameters are varied in turn, one at a time, by adding a fixed proportionality increment to the trajectory baseline values, the model re-run,
#' and the output criterion computed. We refer to these as 'level runs'. The proportionality increment is the same for all of the level runs within a given trajectory, and is drawn at random from a set of fixed
#' levels distributed symetrically around 0 (e.g. -10, -5, +5, +10 percent, i.e. proportions of the trajectory baseline values = 0.9, 0.95, 1.05, 1.10).
#'
#' For each level run, the 'Elementary Effect (EE)' of the given parameter is calculated from the difference between the level run criterion value and the corresponding trajectory baseline criterion value.
#'
#' On completion of all the trajectories, the raw results are (optionally) post-processed to generate the mean and standard deviations of all the EE values for each parameter. EE_mean is an index of the magnitude
#' of the sensitivity, and EE_sd is an index of the extent of interaction with other parameters.
#'
#' During the run the function produces a real-time plot for each trajectory, in which the x-axis represents the sequence of parameters, and the y-axis is the likelihood
#' of the target data. A horizontal red line indicates the likelihood of the parent parameter set, horizontal grey line indicates the likelihood for each trajectory baseline
#' and each level-run likelihood is shown by a symbol. The y-axis range can be changed in real-time by editing the setup file "/Models/Modelname/Modelvariant/Param/control/sensitivity.csv"
#'
#' The outputs from the function are saved as list objects and directed to csv files (provided that the argument csv.output=TRUE) in the "results" folder sepcified in an e2e_read() function call. The outputs are:
#' \itemize{
#'   \item Table of parameter values applied in each run of the model (OAT_parameter_values-*.csv, where * = model.ident as defined by e2e_read())
#'   \item Table of the criterion value and EE value for each trajectory/level run (OAT_results-*.csv)
#'   \item If post-processing is selected, then a table of Mean EE and standard deviation of EE for each parameter, sorted by the absolute value of EE_mean (sorted_parameter_elemental_effects-*.csv)
#' }
#'
#' As mentioned above, the default criterion for assessing the model sensitivity is the likelihood of the observed target data set on the state of the ecosystem given each set of model drivers and parameters. 
#' However, a function argument allows other criteria to be chosen as the basis for the analysis from the list of annually averaged or integrated variables saved in the output objects:
#' \itemize{
#'   \item results$final.year.output$mass_results_wholedomain (whole-domain annual averages of stage variables over the final year of a model run), and 
#'   \item results$final.year.output$annual_flux_results_wholedomain (whole-domian annual integrals of fluxes between state variables over the final year of a model run).
#' }
#' The criterion is chosen by setting a value for the argument outID. The default outID=0 selects the likelihood of the observed target data. Other values in the range 1 to 247 select annualy averaged mass or 
#' annually integrated flux outputs from a list viewable by running the function e2e_get_senscrit().
#' 
#' WARNING - The e2e_run_sens() function will take several days to run to completion on a single processor with even a modest number of iterations. 
#' The total number of model runs required to support the analysis is r*(n+1) where r is the number of trajectories and n is the number of parameters. 
#' The function incorporates all of the physical configuration parameters, fixed and fitted ecology model
#' parameters, the fishing fleet model parameters, and the environmental forcings into the analysis, so n = 450. Each model run needs to be sufficiently long to achieve a stationary state 
#' and as a consequnce a typical runtime will be around 10h per trajectory. The mininimum recommended number of trajectories
#' is 15, so the function can take several days to complete.
#'
#' However, it is possible to spread the load over multiple processor/machines with 
#' arguments in the function allowing for management of this parallelization. Afterwards, the raw results files are combined into a single data set
#' using the e2e_merge_sens_mc() function, and then processed using the function e2e_process_sens_mc(). 
#'
#' A separate function e2e_plot_sens_mc() produces a graphical representation of the EE_mean and EE_sd results.
#'
#' @references For details on the sensitivity analysis method see: Morris, M.D. (1991). Factorial sampling plans for preliminary computational experiments. Technometrics, 33, 161-174.
#' @references For a review of sensitivity analysis methods including the Morris Method see: Wu, J. et al. (2013). Sensitivity analysis of infectious disease models: methods, advances and their application. J R Soc Interface 10: 20121018, 14pp. 
#'
#' @param model R-list object generated by the e2e_read() function which defines the parent model configuration.
#' @param nyears Number of years to run the model in each iteration (default=50).
#' @param n_traj Number of trajectories of parameter sets (default=16).
#' @param trajsd Coefficient of variation used to set the standard deviation for the gaussian distribution from which new paramater values are drawn to create each trajectoy baseline from the initial parent values (default=0.0075).
#' @param n_setoflevels Number of fixed levels of coefficient of variation used to generate the individual parameter values in each level-run. Must be an even number (default=4).
#' @param v_setoflevels Maximum coefficient of variation for the set of levels (default=0.1, i.e. -10 percent to +10 percent).
#' @param coldstart Logical. If TRUE then the run is starting from cold - which means that the first trajectory baseline is the parent configuration as specified in the 'model' list object. If FALSE then signifies that this is a parallel run whuch will later be merged with the 'coldstart=TRUE' run. In this case the first trajectory baseline is a derivative of the parent. Default=TRUE.
#' @param quiet Logical. If TRUE then suppress informational messages at the start of each iteration (default=TRUE).
#' @param postprocess Logical. If TRUE then process the results through to a final sorted list of parameter sensitivies for plotting. If FALSE just produce the raw results. The reason for NOT processing would be if the job has been shared across multiple machines/processors and several raw result files need to be merged before processing. Default=TRUE.
#' @param csv.output Logical. If TRUE then enable writing of csv output files (default=FALSE).
#' @param runtime.plot Logical. If FALSE then disable runtime plotting of the progress of the run - useful for testing (default=TRUE).
#' @param outID Numeric value in the range 0 to 247. Selects the output criterion to be used as the basis for the analysis. Default=0 corresponds to the likelihood of observed target data. Other values obtainable by running e2e_get_senscrit().
#'
#' @return Depends on the settings of arguments 'postprocess' and 'csv.ouptut': If postprocess=TRUE and csv.output=TRUE then outputs are csv files of raw parameter vectors, likelihoods and Elementary Effects for each run, and parameter list sorted by EE_mean, plus the function returns the data of sorted parameters. If postprocess=FALSE and csv.output=FALSE then the function simply returns a dataframe of likelihoods and Elementary Effects for each run.
#'
#' @seealso \code{\link{e2e_get_senscrit}}, \code{\link{e2e_read}}, \code{\link{e2e_merge_sens_mc}}, \code{\link{e2e_process_sens_mc}}, \code{\link{e2e_plot_sens_mc}}
#'
#' @importFrom stats runif 
#' @importFrom grDevices dev.off
#'
#' @export
#'
#' @examples
#' # The examples provided here are illustration of how to set up and run sensitivity
#' # analyses. Even though they are stripped-down minimalist examples, they each still
#' # take about 45 min to run.
#'
#' # --------------------------------------------------------------------------
#'
#' \dontrun{
#' # Load the 2003-2013 version of the North Sea model supplied with the package:
#'     model <- e2e_read("North_Sea", "1970-1999")
#' #
#' # Run the sensitivity analysis process (a quick demonstration):
#' # WARNING - Running a full sensitivity analysis takes days of computer time on a single
#' # machine/processor because it involves a huge number of model runs.
#' # The example below is just a (relatively) quick minimalist demonstration and should NOT
#' # be taken as the basis for any analysis or conclusions. 
#' # Even so, this minimalist demonstration run could take 45 min to complete because it
#' # involves 1353 model runs.
#' # This examples uses the likelihood of observed target data as the criterion for assessing
#' # model sensitivity (i.e. outID is not set and so assumes the default value of 0).
#'     sens_results <- e2e_run_sens(model, nyears=1, n_traj=3, csv.output=FALSE)
#' # View the top few rows of the results dataframe:
#'     head(sens_results)
#' }
#' 
#' # --------------------------------------------------------------------------
#'
#' # To view the list of available criteria for use as the basis for the sensitivity analysis:
#'     e2e_get_senscrit()
#'
#' # --------------------------------------------------------------------------
#'
#' \dontrun{
#' # This examples uses the annual average mass of planktivorous fish as the criterion for
#' assessing model sensitivity (i.e. outID = 24).
#'     sens_results <- e2e_run_sens(model, nyears=1, n_traj=3, csv.output=FALSE, outID=24)
#' # View the top few rows of the results dataframe:
#'     head(sens_results)
#' }
#'
#' # --------------------------------------------------------------------------
#'
#' # This is a dummy example to illustrate a more realistic sensitivity analysis:
#' #    sens_results <- e2e_run_sens(model, nyears=50, n_traj=16, postprocess=TRUE,csv.output=TRUE)
#' # DO NOT launch this configuration unless you are prepared to wait many days for the results
#' 
#' # --------------------------------------------------------------------------
#'
#' \dontrun{
#' # Example of parallelizing the process:
#' # Launch two (or more) runs separately on different processors, with results directed to a 
#' # temporary folder. Set results.path="Yourfolder" to run operataionally. 
#' # Launch batch 1 (on processor 1):
#'     model1 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH1")
#'     sens_results1 <- e2e_run_sens(model1, nyears=1, n_traj=3, coldstart=TRUE, 
#'                    postprocess=FALSE, csv.output=TRUE)
#' # Note that coldstart=TRUE for the first batch only.
#' # Launch batch 2 (on processor 2):
#'     model2 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH2")
#'     sens_results1 <- e2e_run_sens(model2, nyears=1, n_traj=3, coldstart=FALSE, 
#'                                 postprocess=FALSE, csv.output=TRUE)
#' # Note that these two runs return only raw data since postprocess=FALSE
#' 
#' # Then, afterwards, merge the two raw results files with text-tags BATCH1 and BATCH2,
#' # and post process the combined file:
#'     model3 <- e2e_read("North_Sea", "1970-1999", model.ident="COMBINED")
#'     processed_data <- e2e_merge_sens_mc(model3, selection="SENS",
#'                     ident.list<-c("BATCH1","BATCH2"), postprocess=TRUE, csv.output=TRUE)
#' # or...
#'     combined_data <- e2e_merge_sens_mc(model3, selection="SENS",
#'                    ident.list<-c("BATCH1","BATCH2"), postprocess=FALSE, csv.output=TRUE)
#'     processed_data <- e2e_process_sens_mc(model3, selection="SENS",
#'                     use.example=FALSE, csv.output=TRUE)
#' 
#' # Plot a diagram of parameter sensitivities from the combined data
#'     e2e_plot_sens_mc(model3, selection="SENS", use.example=FALSE)
#' }
#'
#' # --------------------------------------------------------------------------
#'
#
# ---------------------------------------------------------------------
# |                                                                   |
# | Authors: Mike Heath, Ian Thurlbeck                                |
# | Department of Mathematics and Statistics                          |
# | University of Strathclyde, Glasgow                                |
# |                                                                   |
# | Date of this version: May 2020                                    |
# |                                                                   |
# ---------------------------------------------------------------------

e2e_run_sens <- function(model, nyears=50, n_traj=16, trajsd=0.0075, n_setoflevels=4, v_setoflevels=0.1, coldstart=TRUE, quiet=TRUE, postprocess=TRUE, csv.output=FALSE, runtime.plot=TRUE,outID=0) {

   oo <- options()
   on.exit(options(oo))

	if(runtime.plot==TRUE){
        start_par = par()$mfrow
	dev.off()
        on.exit(par(mfrow = start_par))
	}

	pkg.env$csv.output <- csv.output	# controls writing of CSV files

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

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

   # Define a function to recompute volume conservation anytime that input flow drivers are jiggled
   vol_conserve <- function(tempforc) {
						#First, these are the flows we actually know
						#		   fdriverso_inflow    (8)
						#		   fdriverd_inflow<-   (9)
						#		   fdriversi_inflow<-  (10)
						#		   fdriversi_outflow<- (13)
						#		   fdriverso_si_flow<- (14)
						#		   fdrivers_upwell<-   (16)

						# The following need to be recomputed every time any o fthese are chnaged in order to maintain volume conservation
						#                  fdriverso_outflow   (11)
						#                  fdriverd_outflow    (12)
						#                  fdriversi_so_flow   (15)

							#Derived outflow from si to so (15)
							#fdriversi_so_flow[,2]<-fdriverso_si_flow[,2] + fdriversi_inflow[,2] - fdriversi_outflow[,2]
							tempforc[[15]][,2]  <- tempforc[[14]][,2] + tempforc[[10]][,2]  - tempforc[[13]][,2]

							#Derived outflow from so (11)
							#fdriverso_outflow[,2]<-fdriverso_inflow[,2] + fdrivers_upwell[,2] + fdriversi_so_flow[,2]  -  fdriverso_si_flow[,2]
							tempforc[[11]][,2]  <- tempforc[[8]][,2] + tempforc[[16]][,2] + tempforc[[15]][,2] - tempforc[[14]][,2]

							#Derived outflow from d (12)
							#fdriverd_outflow[,2]<-fdriverd_inflow[,2] - fdrivers_upwell[,2] 
							tempforc[[12]][,2]  <- tempforc[[9]][,2] - tempforc[[16]][,2] 

		return(tempforc)
   }

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

        # Set up the location of the output metric to be used for the sensitivity analysis
        if(outID<0 || outID > 247) {stop("Error: unknown output metric selection '", outID, "' !\n") }
	if(outID==0) {output_object<-"annual_obj"}
	if(outID>0 & outID<32) {output_object<-"mass_results_wholedomain"}
	if(outID>31) {output_object<-"annual_flux_results_wholedomain"}

	if(outID==0) {rowID<-1}
	if(outID>0 & outID<32) {rowID<-outID}
	if(outID>31) {rowID<-outID-31}

	if(outID==0) {metricname<-"likelihood"}
	if(outID>0 & outID<32) {metricname<-mass_results_descriptions[rowID]}
	if(outID>31) {metricname<-annual_flux_descriptions[rowID]}

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


	#Set the cv level scaling factor to be applied to each parameter
	n_levels<-1		# DONT CHANGE THIS
	#n_setoflevels<-4	# this must be an even number 
	#v_setoflevels<-0.1	#levels are set a +/- this cv arouund zero

	pcvscaleL<-NULL
	for(jj in 1:(n_setoflevels/2)){
		pcvscaleL<-c(pcvscaleL,(1-v_setoflevels)+((jj-1)*v_setoflevels/(n_setoflevels/2)))
	}

	pcvscaleU<-NULL
	for(jj in 1:(n_setoflevels/2)){
		pcvscaleU<-c(pcvscaleU,(1+v_setoflevels)-((jj-1)*v_setoflevels/(n_setoflevels/2)))
	}
	pcvscaleU<-pcvscaleU[seq((length(pcvscaleU)),1,by=-1)]

	pcvscale<-c(pcvscaleL,pcvscaleU)


	#PRODUCES (PARAMETERS+1)*n_traj RUNS - SO BE CAREFUL !

	#End of model setup


	#SO NOW WE NEED TO LOAD THE BEST FIT PARAMETER VECTOR, AND THEN
	#LOOP THROUGH A SET OF RUNS WHERE THE BASELINE VECTOR IS GIVEN A RANDOM JIGGLE EACH TIME
	#SO THIS IS A BIT LIKE THE ANNEALING CODE VERSION WHERE THE ACCEPTED PARAMETERS NEVER CHANGE

	#Vector of the locations in parms which are to be tweaked
	p_to_tweak<-c(seq(1,4,by=1),seq(28,56,by=1),		# physical config parameters excl the area proportions of sediments
		seq(57,123,by=1),seq(286,306,by=1),seq(307,585,by=1))  # FISHING PARAMETERS - omitting the undersize and discard switch parameters loc 124 and 125
                                                                     # and the proportional spatial distributions of discards and offal as these can easily be jiggled
                                                                     # ECOLOGY PARAMETERS - start from end of fleet vector slot through to end of parms
	Nruns<-length(p_to_tweak)+1

	#Set up a vector of parms locations which are constrained to lie within 0 and 1
	p_zero_one<-c(83:123,  # discard and gutting rates
              286:293, # swept area ratios
              294:297, # benthois damage rates
              298,     # offal as prop live weight
              299:306, # plough depths as prop sed layer thickness
	#--------------------------
              480:493, # assimilation efficiencies
              494:507, # excretion rates
              508    , # mt
              509    , # nst
              510    , # dst
              511    , # ndt
              512    , # ddt
              514    , # qs_p2
              515    , # qs_p3
              516    , # msedt
              518    , # nsedt
              520    , # dsedt
              539    , # kelpdebris_det
              540    , # corp_det
              541    , # disc_corp
              542    , # dsink_s
              543    , # dsink_d
              545:550, # migcoef rates
              570:579, # maximum exloitable fractions
              580:583) #fecundities

	#Set up a vector of parms locations where the parameter is expected to be negative
	p_negative<-c(517,519)   # msens,  nsens


        # Vector of drivers to be tweaked - excludes the computed outflows for volume conservation
        # as these are all derived from the input drivers
        f_to_tweak <- c(1:10,13,14,16:53)
	f_additive <- c(4,5,6) # temperatures

	whichforclogs<-c(2,3)

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

	build <- build_model(model ,nyears)
	parms <- elt(build, "model.parameters")
	parmname<-(names(parms))[p_to_tweak]

	forc <- elt(build, "forcings")

	Nforcstot<-length(names(forc))
	Nforcs<-length(f_to_tweak)
	forcid<-(seq(1001,1000+Nforcstot,by=1))[f_to_tweak]
	forcname<-(names(forc))[f_to_tweak]

	#add Nforcs to Nruns
	Nruns<-Nruns+Nforcs

	def_parms <- parms
	def_forc <- forc


	#ITERATE THE TRAJECTORY PROCESS.....


	#Set up an array to hold the results which will be the selected metric value from the the model
	OAT_results<-array(NA,dim=c(n_levels*Nruns*n_traj,6))
	colnames(OAT_results)<-c("parameterid","trajectoryid","levelid","delta_p",metricname,"EE")
	rownames(OAT_results)<-rep(rep(c("baseline",forcname,names(def_parms[p_to_tweak])),n_levels),n_traj)

	OAT_results[,1]<-rep(c(0,forcid,p_to_tweak),n_traj*n_levels)

	traji<-rep(1,Nruns*n_levels)
	traj<-NULL
	for(iiz in 1:(n_traj)) {
		traj<-c(traj,(traji+iiz-1))
	}

	OAT_results[,2]<-traj

	levid<-NULL
	lev<-NULL
	for(iiz in 1:(n_traj)) {
		level_to_use<-floor(runif(1,1,(n_setoflevels + 1 -(1e-10))))
		levid<-c(levid,rep(level_to_use,Nruns))
		lev<-c(lev,rep((pcvscale[level_to_use]-1),Nruns))
	}	# NOTE (p-1) here

	OAT_results[,3]<-levid
	OAT_results[,4]<-lev
	OAT_results[which(OAT_results[,1]==0),3]<-0
	OAT_results[which(OAT_results[,1]==0),4]<-0

	#Set up an array to hold the vectors of parameter values
	OAT_parmvalues <- array(NA,dim=c((Nruns)*n_levels*n_traj,(Nruns-Nforcs)-1))
	colnames(OAT_parmvalues)<-names(def_parms[p_to_tweak])


        trigger<-1
        if(coldstart==FALSE) trigger<-0

	tot_exp_runs <- n_traj * n_levels * Nruns   # expected total number of runs		

	for (kkk in 1:n_traj){

		traj_parms<-def_parms
		traj_forc<-def_forc

		if(kkk>trigger){
			# FIRST RUN FROM A COLD START USES THE MAX LIKELIHOOD AS THE BASELINE SO SET kkk>1
			# IF THIS IS A HOT START THEN CHANGE THIS TO kkk>0

			#Jiggle all the forcings to create a new trajectory

			for(ffff in f_to_tweak){
	
				if(length( which(whichforclogs==ffff) ) == 0){
					if(length( which(f_additive==ffff) ) == 0){       # multiplicative scaling
					deltaforc<-rnorm(1,1,trajsd)
					traj_forc[[ffff]][,2]<-def_forc[[ffff]][,2]*deltaforc
					}
					if(length( which(f_additive==ffff) ) > 0){        # additive scaling
					deltaforc<-rnorm(1,1,trajsd)
					traj_forc[[ffff]][,2]<-def_forc[[ffff]][,2] + mean(def_forc[[ffff]][,2]) * (deltaforc-1)
					}
				}

				if(length( which(whichforclogs==ffff) ) > 0){
					deltaforc<-rnorm(1,1,trajsd)
					traj_forc[[ffff]][,2]<-(def_forc[[ffff]][,2]) + log(deltaforc)
				}

			}


			#Now we have to be a bit careful and redo the volume conservation calculations
			#First, these are the flows we actually know
			#		   fdriverso_inflow    (8)
			#		   fdriverd_inflow<-   (9)
			#		   fdriversi_inflow<-  (10)
			#		   fdriversi_outflow<- (13)
			#		   fdriverso_si_flow<- (14)
			#		   fdrivers_upwell<-   (16)
                        traj_forc <- vol_conserve(traj_forc)


			#Jiggle all the parameters (inlcuding the fixed ones) to create a a new trajectory

			for(aaaaa in p_to_tweak){

				if(length( which(p_negative==aaaaa) ) == 0){
					traj_parms[aaaaa]<-max(0,rnorm(1,def_parms[aaaaa],def_parms[aaaaa]*trajsd))
				}

				if(length( which(p_negative==aaaaa) ) > 0){
					traj_parms[aaaaa]<-max(0,rnorm(1,(-1*def_parms[aaaaa]),(-1*def_parms[aaaaa])*trajsd))
					traj_parms[aaaaa]<-(-1*traj_parms[aaaaa])
				}

				#Then any that are not allowed to be greater than 1 are trimmed back to 1
				if(length( which(p_zero_one==aaaaa) ) >0){
					if(traj_parms[aaaaa]>1) traj_parms[aaaaa]<-1
				}

			}

		}

		#===========================================


		#So, here we go....
		for(qqqq in 1:n_levels){	# for each level

			for(jjjj in 1:Nruns){	# for each parameter

			annealing_control_data<-read_SD_sensitivity(model.path)

			if(jjjj==1){
			axmin <- elt(annealing_control_data, "axmin")
			axmax <- elt(annealing_control_data, "axmax")
			axmin_p<-axmin
			axmax_p<-axmax
			newglims<-1
			}
			if(jjjj>1){

                        baseline_traj_value <- OAT_results[1,5]

			newglims<-0
			axmin <- elt(annealing_control_data, "axmin")
				if(baseline_traj_value < axmin) {
				message("Warning: Baseline < axis min in sensitivity.csv control file - axmin reset to zero")
				axmin<-0
				}
				if((axmin==axmin_p)==FALSE){
					newglims<-1
				}
			axmax <- elt(annealing_control_data, "axmax")
				if(baseline_traj_value > axmax) {
				message("Warning: Baseline > axis max in sensitivity.csv control file - axmax reset to baseline*1.5")
				axmax<-baseline_traj_value*1.5
				}
				if((axmax==axmax_p)==FALSE){
					newglims<-1
				}
			axmin_p<-axmin
			axmax_p<-axmax
			}

			if( axmax<=axmin ) {
			message("Warning: Plot axis max <= axis min in sensitivity.csv control file - reset to defaults")
			axmax<-1
			axmin<-0
			newglims<-1
			}

				rtof<-(Nruns*n_levels*(kkk-1)) + (Nruns*(qqqq-1)) + jjjj
				rtst<-(Nruns*n_levels*(kkk-1)) + (Nruns*(qqqq-1)) + 1
				rten<-(Nruns*n_levels*(kkk-1)) + (Nruns*(qqqq-1)) + Nruns

				forc<-traj_forc

				parms<-traj_parms

				if(jjjj>1) {

					if(jjjj<(Nforcs+2)){

				forc_to_jig <- f_to_tweak[jjjj-1]

				if(length( which(whichforclogs==forc_to_jig) ) == 0){
					if(length( which(f_additive==forc_to_jig) ) == 0){       # multiplicative scaling
					deltaforc<-pcvscale[OAT_results[rtof,3]]
					forc[[forc_to_jig]][,2]<-traj_forc[[forc_to_jig]][,2]*deltaforc
					}
					if(length( which(f_additive==forc_to_jig) ) > 0){        # additive scaling
					deltaforc<-pcvscale[OAT_results[rtof,3]]
					forc[[forc_to_jig]][,2]<-traj_forc[[forc_to_jig]][,2] + mean(traj_forc[[forc_to_jig]][,2]) * (deltaforc-1)
					}
				}

				if(length( which(whichforclogs==forc_to_jig) ) > 0){
					forc[[forc_to_jig]][,2]<-(traj_forc[[forc_to_jig]][,2]) + log(pcvscale[OAT_results[rtof,3]])
				}

					# ........................................................................

						#Now we have to be a bit careful and redo the volume conservation calculations

						#if jjjj is in the range 8 to 16 then we need to redo the volume conservation calculation
						#These are the flows we actually know
						#		   fdriverso_inflow    (8)
						#		   fdriverd_inflow<-   (9)
						#		   fdriversi_inflow<-  (10)
						#		   fdriversi_outflow<- (13)
						#		   fdriverso_si_flow<- (14)
						#		   fdrivers_upwell<-   (16)
						# The following need to be recomputed every time any o fthese are chnaged in order to maintain volume conservation
						#                  fdriverso_outflow   (11)
						#                  fdriverd_outflow    (12)
						#                  fdriversi_so_flow   (15)

						if(length(which( (c(8,9,10,13,14,16))==(forc_to_jig))  ) > 0) forc <- vol_conserve(forc)
   					}


   					if(jjjj>(Nforcs+1)) {
						# ZZ we are writing the parms used in the model run here:
   						parms[p_to_tweak[jjjj-Nforcs-1]] <- traj_parms[p_to_tweak[jjjj-Nforcs-1]] * pcvscale[OAT_results[rtof,3]]
   						if(length(     which(p_zero_one==parms[p_to_tweak[jjjj-Nforcs-1]])) >0 ) {
      							if(parms[p_to_tweak[jjjj-Nforcs-1]] > 1) parms[p_to_tweak[jjjj-Nforcs-1]]<-1
   						}
	   				}

				}

				needtorun<-1
				if(jjjj>(Nforcs+1)){
					if(def_parms[p_to_tweak[jjjj-Nforcs-1]]==0){
     						needtorun<-0  # if the default parameter value is zero no need to run the model
					}
				}

				if(needtorun==1){
					# need to put modified forcings back into build:
					build$forcings <- forc
					build$model.parameters <- parms


########################
# Select which model run function to use and then extract the approriate model output metric

	if(outID==0)   {
		fit <- StrathE2E_fit(model, build, quiet)
		Partial_chi <- elt(fit, "partial_chi")
		metric_value <- Partial_chi[nrow(Partial_chi),]
	}

	if(outID>0)   {
	output <- StrathE2E.run(build)
	aggregates			<- aggregate_model_output(model, output)
	annual.results.wholedomain	<- derive_annual_results_wholedomain(model, build, output, aggregates)
	if(outID>0 & outID<32)      { metric_value <- annual.results.wholedomain$mass_results[rowID,1] }
	if(outID>31)                { metric_value <- annual.results.wholedomain$annual_flux_results[rowID,1] }
	}

#########################


				}

				OAT_parmvalues[rtof,]<-parms[p_to_tweak]

#INSERT THE REQUIRED SENSITIVITY OBJECT INTO OAT_results[rtof,5]

				if(needtorun==1) OAT_results[rtof,5] <- metric_value
				if(needtorun==0) OAT_results[rtof,5] <- OAT_results[rtst,5]		# if default parm value is zero implant the baseline metric value

				if(jjjj==1) OAT_results[rtof,6]<-0
				if(jjjj>1) OAT_results[rtof,6]<- ( (OAT_results[rtof,5]-OAT_results[rtst,5])/ sqrt((OAT_results[rtof,4]^2)) )

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

		#Plot or update the time series of trajectory baseline and individual metric values so far....

				if(runtime.plot==TRUE){
				if(jjjj>1) {

					if(jjjj==2 | newglims==1){
					par(mfrow=c(1,1))
					pltxset<-c( (forcid) , p_to_tweak )
					pltxtyp<-c( 	rep("black",(length(forcid))),
							rep("red",(33-1+1)),
							rep("blue",(121-34+1)),
							rep("green",(400-122+1)) )
					pltxlabels<-c(	"Environmental drivers",
							"Physical setup parameters",
							"Fishing related parameters",
							"Ecology parameters",
							"Initial trajectory baseline",
							"Current trajectory baseline")
					plot(c(1,(Nruns-1)) , c(OAT_results[1,5],OAT_results[1,5]) , type="l", col="red", ylim=c(axmin,axmax), xlim=c(1,(Nruns-1)),
						xlab="Parameters",ylab=metricname,
						main=paste("Trajectory ",kkk,sep="") )
					points( c(1,(Nruns-1)) , c(OAT_results[rtst,5],OAT_results[rtst,5]),type="l" , col="grey")
				        points(seq(1:(jjjj-1)),OAT_results[(rtst+1):rtof,5],type="p",pch=20, col=pltxtyp[1:(jjjj-1)])
			                legend("topleft",pltxlabels,bg="transparent",
						col=c("black","red","blue","green","red","grey"),
						lty=c(0,0,0,0,1,1),
						lwd=c(0,0,0,0,1,1),
						pch=c(20,20,20,20,NA,NA),
						pt.cex=c(1,1,1,1,0,0),
						cex=c(0.7,0.7,0.7,0.7,0.7,0.7))
					}

					if(newglims==0){ 
				        points((jjjj-1),OAT_results[rtof,5],type="p",pch=20, col=pltxtyp[(jjjj-1)])
					}

				}
				}

#----------------------------------------------------------------------------------------------------------
				
				pc_complete<-  (rtof*100/tot_exp_runs)
				message(sprintf("%.3f",pc_complete),"% complete; current trajectory ",kkk," of ",n_traj,"; parameter ",jjjj," of ",Nruns,"           ",appendLF=TRUE)

				#-------------------
				#Transfer the results into a dataframe for output (so number of elements in the csv header row = number of columns)
				# - makes it easier to read the data in again
				results_df_out <- as.data.frame(rep(0,nrow(OAT_results)))
				names(results_df_out) <-"parametername"
				results_df_out$parameterid<-OAT_results[,1]
				results_df_out$trajectoryid<-OAT_results[,2]
				results_df_out$levelid<-OAT_results[,3]
				results_df_out$delta_p<-OAT_results[,4]
				results_df_out$outputmetric<-OAT_results[,5]  # renamed later to the dynamic value selected by function argument outID
				results_df_out$EE<-OAT_results[,6]
				results_df_out$parametername <- rownames(OAT_results)

	                        names(results_df_out)<-c("parametername","parameterid","trajectoryid","levelid","delta_p",metricname,"EE")

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

				filename <- csvname(resultsdir, "OAT_results", model.ident)
				writecsv(results_df_out, filename, row.names=FALSE)

				filename <- csvname(resultsdir, "OAT_parameter_values", model.ident)
				writecsv(OAT_parmvalues, filename, row.names=FALSE)

			} # end inner Runs jjjj loop
		} # end mid levels qqqq loop
	} # end outer trajectories kkk loop

	message("")

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

	#PROCESS THE FINAL RESULTS TO GET THE SORTED PARAMETER SENSITIVITY FILE IF REQUIRED
	if(postprocess==TRUE){
	message("Processing results... ", appendLF=FALSE)
		output_data <- process_sensitivity_analysis_results(model, results_df_out)
	message("Completed")
        }

	if(postprocess==FALSE){
		output_data <- results_df_out
        }

	output_data

}