R/optimize_activity_mult_eco.R

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

#
# optimize_activity_mult_eco.R
#
#' Optimise StrathE2E fishing gear activity multipliers to maximum the likelihood of observed ecosystem target data.
#'
#' Launches a simulated annealing process to find the set of fishing fleet model gear activity multipliers
#' producing the maximum likelihood of observed target data on the state of the ecosystem,
#' given specified environmental driving data and ecology model parameters, and effort-harvest ratio scaling parameters.
#' Note that there is a related function optimize_activity_mult_hr() which uses only the fishing fleet model to
#' find the activity multipliers which maximise the likelihood of a set of target harvest ratios (where these are known), given values of the
#' scaling values linking effort to harvest ratio.
#'
#' 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)
#'
#' Model configuration and initial values of the ecology model parameters need to be
#' assembled by a prior call of the e2e_read() function.
#' 
#' NOTE that the models.path argument in the e2e_read() function call needs to point to a user workspace folder, not the default
#' North Sea model provided with the package. This is because the annealing function needs write-access to the model /Param folder,
#' but the /extdata/Models folder in the package installation is read-only.
#' To use the annealing function on the North Sea model, use the e2e_copy() function to make a copy of the
#' North Sea model in the user workspace.
#'
#' The coefficient of variation for jiggling the activity multipliers can be varied in real-time
#' during the run by editing the file "optimize_fishing.csv" in the folder
#' /Param/control/ of the model version. Suggested vaues for the CV are in the range 0.1 to 0.01
#'
#' The function produces a real-time graphical summary of the progress of the fitting procedure, displaying
#' the likelihoods of the proposed and accepted parameter sets at each iteration.
#' Y-axis (likelihood of the target data) range of the real time plot can be varied during the run by 
#' editing the setup file "optimize_fishing.csv"
#'
#' At the end of the procedure a new version of the gear activity multipliers file is exported
#' to the folder /Param of the model version, with a user defined identifier specified by the model.ident argument
#' in the e2e_read() function. The histories of proposed and accepted parameter combinations
#' are saved as csv files in the results folder.
#'
#' To preserve the new activity multipliers and incorporate them into the fishing fleet model parameterisation
#' the multiplier values need to be applied to the activity rates specified in the data input file /Param/fishing_activity_*.csv.
#' Manually update the values in fishing_activity_*.csv, by multiplying the existing values by the new multipliers emerging from the annealing process.
#'
#' If the edited file fishing_activity_*.csv is saved with a new identifier (*) then in order to use it in a subsequent run of the
#' StrathE2E model (using the e2e_run() function) it will be necessary to edit the MODEL_SETUP.csv file in the relevant /Models/variant folder to point to the new file.
#'
#' @param model R-list object generated by the e2e_read() function which defined the model configuration.
#' @param nyears Number of years to run the model in each iteration (default=40).
#' @param n_iter Number of iterations of the model (default=500).
#' @param start_temperature Initial value of the simulated annealing temperature parameter (default=1). Suggested values in the range 0.0005 - 5. Higher values increase the probability of rejecting parameter combinations producing an improvement in likelihood.
#' @param cooling Rate at which the simulated annealing temperature declines with iterations (default=0.975). Suggested values in the range 0.9 - 0.985
#' @param quiet Logical. If TRUE then suppress informational messages at the start of each iteration, (default=TRUE).
#' @param csv.output Logical. If FALSE then disables writing of csv output files - useful for testing, (default=TRUE)
#' @param runtime.plot Logical. If FALSE then disables runtime plotting of the progress of the run - useful for testing, (default=TRUE).
#'
#' @return A list object containing the histories of proposed and accepted parameter values and the final accepted parameter values. Optionally (by default), CSV files of the histories and the final accepted parameter values. The latter are returned to the model parameter folder in a format to be read back into the model.
#'
#' @importFrom stats runif
#' @importFrom grDevices dev.off
#'
#' @noRd
#
# ------------------------------------------------------------------------------

optimize_activity_mult_eco <- function(model, nyears=40, n_iter=500, start_temperature=1, cooling=0.975, quiet=TRUE, csv.output=TRUE, runtime.plot=TRUE) {

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

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

	data				<- elt(model, "data")
	fleet.model			<- elt(data, "fleet.model")
	fitted.parameters		<- elt(data, "fitted.parameters")
	gear_mult			<- elt(fleet.model, "gear_mult")
	HRscale_vector_multiplier	<- elt(fleet.model, "HRscale_vector_multiplier")

	if (read.only & csv.output==TRUE) {
		message("Warning: cannot write fitted parameters back to the model input folders - model is read-only")
		message("Warning: to fix this, make a copy of the model using e2e_copy() into your own workspace.")
		stop("Model is not writable!")
	}

	print(date())


	if (length(which(HRscale_vector_multiplier<1 | HRscale_vector_multiplier>1)) > 0) {
		print("**************************************************************************")
		print("WARNING - one or more baseline Harvest Ratio scaling values differs from 1")
		print("**************************************************************************")
	}

	#Set up a dataframe to hold the proposed and accepted ACT_scale_vector

	ACTmult_proposal<-data.frame(G1=0,G2=0,G3=0,G4=0,G5=0,G6=0,G7=0,G8=0,G9=0,G10=0,G11=0,G12=0,lik=0)
	ACTmult_proposal[1,]<-c(gear_mult,1e-60)
	ACTmult_accepted<-ACTmult_proposal

	n_acceptances <- 0	# Counter for the number of parameter acceptances which have occurred
	lastbest <- 1e-60
	
	temperature <- start_temperature

	#ITERATE THE ANNEALING PROCESS.....

                for (kkk in 1:n_iter){

		annealing_control_data <- read_SD_fishing(model.path)

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

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

		deltaH <- elt(annealing_control_data, "deltaH")		# cv for jiggling the fishing fleet parameters
		#Values for the cv are held in the "optimize_fishing" file.
		#THESE CAN BE ALTERED DURING A RUN BY MANUALLY EDITING THE VALUES IN THE CONTROL FILE
		# Suggested values are:
		#    deltaH<-0.10  # high value for the first cycle through all this
		#    deltaH<-0.05  # medium value is strating from something half sensible
		#    deltaH<-0.02  # small vakue for the second cycle through all this
		#    deltaH<-0.01  # very small value if strating from something close

		if(kkk>1){
			for(www in 1:12){
				gear_mult[www] <- max(0,rnorm(1,ACTmult_accepted[(kkk-1),www],deltaH*ACTmult_accepted[(kkk-1),www]))
				#Might be useful in some cases: force the MF HR scaler to be same as PF HR scaler
				#if(www==3) HRscale_vector_multiplier[3]<-HRscale_vector_multiplier[1]
			}
			ACTmult_proposal[kkk,]<-c(gear_mult,0)
			ACTmult_accepted[kkk,]<-ACTmult_accepted[(kkk-1),]

			#Now simply drop the new fleet vector into the existing parms vector
			model$data$fleet.model$gear_mult <- gear_mult
		}

		# showall("gear_mult", gear_mult)

		# Build and run the model:
		build <- build_model(model, nyears)
		fit <- StrathE2E_fit(model, build, quiet)
		annual_obj <- elt(fit, "annual_obj")

		# message("Err=", annual_obj)

		ACTmult_proposal[kkk,13]<-annual_obj

		temperature<-temperature*cooling

		#-------------------------------------------------------------------
		#Now the Metropolis algorithm.....

		lik_ratio<-exp(((log(annual_obj)) - log(lastbest))/temperature)
		rand<-runif(1,0,1)

		new_accepted<-" accepted: NO"
		if(lik_ratio>rand){
			n_acceptances<-n_acceptances+1
			ACTmult_accepted[kkk,]<-ACTmult_proposal[(kkk),]
			lastbest<-annual_obj

			new_accepted<-" accepted: YES"
		}

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

		filename = csvname(resultsdir, "annealing_ACTmult_proposalhistory", identifier)
		writecsv(ACTmult_proposal, filename, row.names=FALSE)

		filename = csvname(resultsdir, "annealing_ACTmult_acceptedhistory", identifier)
		writecsv(ACTmult_accepted, filename, row.names=FALSE)

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

		message("Iteration: ",kkk,"; proposal likelihood: ",sprintf("%.7f",ACTmult_proposal[kkk,13]),"; ",new_accepted)

		#print(paste("Itn:",kkk," ",ACTmult_proposal[kkk,13]," ",new_accepted,sep=""))
		#print(ACTmult_proposal[kkk,1:12])

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

		#Plot or update the time series of proposal and acepted likelihoods so far....
		if(runtime.plot==TRUE){
		if(kkk>1){

			if(kkk==2 | newglims==1){
			par(mfrow=c(1,1))
			plot(seq(1,nrow(ACTmult_proposal)),ACTmult_proposal[,13],ylim=c(axmin,axmax),xlim=c(1,n_iter),xlab="Iterations",ylab="Target data likelihood",type="l",col="grey")
			points(seq(1,nrow(ACTmult_accepted)),ACTmult_accepted[,13],type="l",col="black",lwd=3)
                	legend("topleft",c("accepted","proposed"),bg="transparent",col=c("black","grey"),lty=c(1,1),lwd=c(3,1),pt.cex=c(1,1))
			} 

			if(newglims==0){ 
			stp<-c((kkk-1),kkk)
			points(stp,ACTmult_proposal[stp,13],type="l",col="grey")
			points(stp,ACTmult_accepted[stp,13],type="l",col="black",lwd=3)
			}

		}
		}

	}

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

	new_parameters<-extract_ACTmult_to_parmsfolder(model, ACTmult_accepted,csv.output=csv.output)

	#Assemble a list object containing all the data generated by the annealing process
	output_data<-list(parameter_proposal_history = ACTmult_proposal,
			  parameter_accepted_history = ACTmult_accepted,
			  new_parameter_data = new_parameters)
	output_data

}