Optimize StrathE2EPolar ecology model parameters to maximise the likelihood of observed ecosystem target data.
e2ep_optimize_eco.RdLaunches a StrathE2EPolar simulated annealing process to find the set of ecology model parameters producing the maximum likelihood of observed target data on the state of the ecosystem, given specified environmental driving data and fishing fleet parameters.
Usage
e2ep_optimize_eco(
model,
nyears = 40,
n_iter = 500,
start_temperature = 1,
cooling = 0.975,
toppredlock = TRUE,
quiet = TRUE,
csv.output = FALSE,
runtime.plot = TRUE
)Arguments
- model
R-list object generated by the e2ep_read() function which defined the model configuration.
- nyears
Number of years to run the model in each iteration (default=40).
- n_iter
Number of iterations of the model (default=500).
- 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.
- cooling
Rate at which the simulated annealing temperature declines with iterations (default=0.975). Suggested values in the range 0.9 - 0.985
- toppredlock
Logical. If TRUE then locks-down the uptake parameters of the birds pinnipeds and cetaceans as these are hard to fit alongside the other parameters (default=TRUE).
- quiet
Logical. If TRUE then suppress informational messages at the start of each iteration (default=TRUE).
- csv.output
Logical. If TRUE then enable writing of csv output files (default=FALSE).
- runtime.plot
Logical. If FALSE then disable runtime plotting of the progress of the run - useful for testing (default=TRUE)
Value
A list object containing the histories of proposed and accepted parameters 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.
Details
Simulated annealing is is a probabilistic technique for approximating the global optimum of a given function. As implemented here the process 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 e2ep_read() function.
NOTE that the models.path argument in the e2ep_read() function call needs to point to a user workspace folder, not the default Barents 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 Barents Sea model, use the e2ep_copy() function to make a copy of the Barents Sea model in the user workspace.
The observational data to which the ecology parameters are optimized are loaded from the folder Modelname/Variantname/Target/annual_observed_*.csv as part of a e2ep_read() function call and are built into the R-list object generated by e2ep_read(). Column 3 of annual_observed_* (header: "Use1_0") is a flag to set whether any given row is used in calculating the likelihood of the observed data given the model setup and parameters. Un-used rows of data are omitted from calculations.
The coefficients of variation for jiggling the ecology parameter can be varied in real-time during the run by editing the file "optimize_ecology.csv" in the folder /Param/control/ of the model version.
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_ecology.csv"
At the end of the procedure, provided that csv.output=TRUE, new versions of the three ecology model 'fitted_parameters..' files are exported to the folder /Param of the model version, with a user defined identifier specified by the model.ident argument in the e2ep_read() function. These data are also saved in the list object returned by the function.
In order to use the new fitted parameter values in a subsequent run of the StrathE2E model (using the e2ep_run() function) it will be necessary to edit the MODEL_SETUP.csv file in the relevant /Models/variant folder to point to the new files.
Also at the end of the procedure the histories of proposed and accepted ecology model parameter values and corresponding likleihoods from each iteration of the procedure are saved as CSV files in the results folder (provided that the argument csv.output=TRUE), and in a list object which is returned by the function. The two csv files generated by the procedure have names: annealing_par_proposalhistory-*, annealing_par_acceptedhistory-*, where * denotes the value of model.ident defined in the preceding e2ep_read() function call. The returned list object contains three dataframes: parameter_proposal_history, parameter_accepted_history, new_parameter_data (a list of three). The proposal and accepted histories can be further analysed with the function e2ep_plot_opt_diagnostics() to assess the performance of the optimization process.
Examples
# \donttest{
# Load the 2011-2019 version of the Barents Sea model supplied with the package and generate a
# quick test data object with only 8 itereations and running the model for only 3 years.
# Also, the final parameter values are not saved back to the model Param folder.
# More realistic would be at least 500 iterations and running for 50 years.
# Even so this example will take a few minutes to run:
model<-e2ep_read(model.name="Barents_Sea",
model.variant="2011-2019",
model.ident="test")
#> Current working directory is...
#> 'C:/Users/jackl/OneDrive - University of Strathclyde/Documents/StrathE2E/strath-e-2-e-polar-webdev/docs/reference'
#> No 'results.path' specified so any csv data requested
#> will be directed to/from the temporary directory...
#> 'C:\Users\jackl\AppData\Local\Temp\RtmpSgdWsc'
#>
#> Model setup and parameters gathered from ...
#> StrathE2E2 package folder
#> Model results will be directed to/from ...
#> 'C:\Users\jackl\AppData\Local\Temp\RtmpSgdWsc/Barents_Sea/2011-2019/'
# This model is already optimized to the observed ecosystem data supplied with the package
# so to illustrate the performance of the process we perturb the temperature driving to knock
# the model away from its maximum likelihood state relative to the target data:
# add 3 degC to upper layer offshore temperatures:
model$data$physics.drivers$so_temp <- model$data$physics.drivers$so_temp+3
# add 3 degC to inshore temperatures:
model$data$physics.drivers$si_temp <- model$data$physics.drivers$si_temp+3
# add 3 degC to lower layer offshore temperatures:
model$data$physics.drivers$d_temp <- model$data$physics.drivers$d_temp+3
test_run <- e2ep_optimize_eco(model, nyears=3, n_iter=8, start_temperature=0.4,
csv.output=FALSE)
#> [1] "Wed Dec 14 21:10:01 2022"
#> Iteration: 1; proposal likelihood: 0.4910439; accepted: YES
#> Iteration: 2; proposal likelihood: 0.4905202; accepted: YES
#> Iteration: 3; proposal likelihood: 0.4907998; accepted: YES
#> Iteration: 4; proposal likelihood: 0.4914530; accepted: YES
#> Iteration: 5; proposal likelihood: 0.4903976; accepted: YES
#> Iteration: 6; proposal likelihood: 0.4897374; accepted: YES
#> Iteration: 7; proposal likelihood: 0.4905885; accepted: YES
#> Iteration: 8; proposal likelihood: 0.4899016; accepted: YES
# View the structure of the returned list:
str(test_run,max.level=1)
#> List of 3
#> $ parameter_proposal_history:'data.frame': 8 obs. of 192 variables:
#> $ parameter_accepted_history:'data.frame': 8 obs. of 192 variables:
#> $ new_parameter_data :List of 3
# View the structure of the returned list element containing parameter objects:
str(test_run$new_parameter_data,max.level=1)
#> List of 3
#> $ new_preference_matrix :'data.frame': 26 obs. of 18 variables:
#> $ new_uptake_mort_rate_parameters:'data.frame': 19 obs. of 9 variables:
#> $ new_microbiology_parameters :'data.frame': 26 obs. of 2 variables:
# View the new preference matrix:
test_run$new_parameter_data$new_preference_matrix
#> kelp icealg phyt omnivzoo carnzoo fishplar
#> ammonia 0.2128508 0.3150904 0.1989438 NA NA NA
#> nitrate 0.7871492 0.6849096 0.8010562 NA NA NA
#> suspdet NA NA NA 0.52234959 NA NA
#> seddet NA NA NA NA NA NA
#> kelpdebris NA NA NA NA NA NA
#> icedet NA NA NA 0.08314314 NA NA
#> corpses NA NA NA NA NA NA
#> discards NA NA NA NA NA NA
#> kelp NA NA NA NA NA NA
#> icealg NA NA NA 0.13359398 NA NA
#> phyt NA NA NA 0.22493457 NA NA
#> omnivzoo NA NA NA NA 0.59417631 0.1645996
#> carnzoo NA NA NA NA NA NA
#> fishplar NA NA NA NA 0.03231941 NA
#> fishdlar NA NA NA NA 0.04046481 NA
#> fishp NA NA NA NA NA NA
#> fishm NA NA NA NA NA NA
#> fishd NA NA NA NA NA NA
#> benthslar NA NA NA 0.01667362 0.23135953 0.5920421
#> benthclar NA NA NA 0.01930510 0.10167993 0.2433583
#> benths NA NA NA NA NA NA
#> benthc NA NA NA NA NA NA
#> bird NA NA NA NA NA NA
#> seal NA NA NA NA NA NA
#> ceta NA NA NA NA NA NA
#> bear NA NA NA NA NA NA
#> fishdlar fishp fishm fishd benthslar benthclar
#> ammonia NA NA NA NA NA NA
#> nitrate NA NA NA NA NA NA
#> suspdet NA NA NA NA 0.98037027 0.2730892
#> seddet NA NA NA NA NA NA
#> kelpdebris NA NA NA NA NA NA
#> icedet NA NA NA NA NA NA
#> corpses NA NA NA 0.007913448 NA NA
#> discards NA NA NA 0.100974259 NA NA
#> kelp NA NA NA NA NA NA
#> icealg NA NA NA NA NA NA
#> phyt NA NA NA NA 0.01962973 0.7269108
#> omnivzoo 0.957975799 0.36806503 0.33556847 NA NA NA
#> carnzoo NA 0.16307892 0.01530567 0.080655042 NA NA
#> fishplar NA 0.12180992 0.06191863 0.149790611 NA NA
#> fishdlar NA 0.15616008 0.08083372 0.096282816 NA NA
#> fishp NA NA NA 0.108528930 NA NA
#> fishm NA NA NA 0.029785019 NA NA
#> fishd NA NA NA 0.006244924 NA NA
#> benthslar 0.034677025 0.09561029 0.27266436 NA NA NA
#> benthclar 0.007347176 0.09527575 0.23370916 NA NA NA
#> benths NA NA NA 0.349951331 NA NA
#> benthc NA NA NA 0.069873619 NA NA
#> bird NA NA NA NA NA NA
#> seal NA NA NA NA NA NA
#> ceta NA NA NA NA NA NA
#> bear NA NA NA NA NA NA
#> benths benthc bird seal ceta bear
#> ammonia NA NA NA NA NA NA
#> nitrate NA NA NA NA NA NA
#> suspdet 0.4090218 NA NA NA NA NA
#> seddet 0.3631364 NA NA NA NA NA
#> kelpdebris NA 0.004367581 NA NA NA NA
#> icedet NA NA NA NA NA NA
#> corpses NA 0.486021267 0.08884347 0.03431177 NA 0.10509925
#> discards NA NA 0.26052450 0.05887910 0.148077608 NA
#> kelp NA 0.009878544 NA NA NA NA
#> icealg NA NA NA NA NA NA
#> phyt 0.2278418 NA NA NA NA NA
#> omnivzoo NA NA NA NA 0.011518176 NA
#> carnzoo NA NA 0.04757628 0.09063836 0.022248727 NA
#> fishplar NA NA NA NA NA NA
#> fishdlar NA NA NA NA NA NA
#> fishp NA NA 0.49229286 0.69922205 0.379060225 NA
#> fishm NA NA 0.10489082 0.04093234 0.397393592 NA
#> fishd NA NA 0.00587208 0.03836033 0.036779705 NA
#> benthslar NA NA NA NA NA NA
#> benthclar NA NA NA NA NA NA
#> benths NA 0.499732607 0.00000000 0.01781157 0.000000000 NA
#> benthc NA NA 0.00000000 0.01984447 0.000000000 NA
#> bird NA NA NA 0.00000000 0.000000000 0.09120423
#> seal NA NA NA NA 0.004921967 0.78269802
#> ceta NA NA NA NA NA 0.02099850
#> bear NA NA NA NA NA NA
# }
# --------------------------------------------------------------------------
# This is a dummy example to illustrate a realistic run in which optimised
# parameters are written back to the model Param folder. To try it out substitute
# your own relative folder path in place of \Folder in the e2ep_copy() function...
# WARNING - this will take about 26 hours to run...
# Copy the 2011-2019 version of the Barents Sea model supplied with the package into a
# user workspace relative to the current working directory (../Folder):
# e2ep_copy("Barents_Sea", "2011-2019",
# dest.path="Folder")
# Load the copied version of the Barents Sea/2011-2019 model from the user workspace
# and assign a path for results data:
# (REPLACE "Folder/Models" and "Folder/results" with your own paths before running)
# model<-e2ep_read(model.name="Barents_Sea",
# model.variant="2011-2019",
# models.path="Folder/Models",
# results.path="Folder/results",
# model.ident="fittingrun")
# Launch the fitting process
# fitting_data <- e2ep_optimize_eco(model, nyears=50, n_iter=500, start_temperature=1,
# csv.output=TRUE)
# --------------------------------------------------------------------------