This notebook tracks progress on the development of GMSE software R package, for game-theoretic management strategy evaluation, and related issues surrounding the development and application of game theory for addressing questions of biodiversity and food security.

# Contents:

Towards a Game-theoretic Management Strategy Evaluation (G-MSE)

General model development

General software development

Game-theory modelling (game.c; green box above)

Game-theory and modelling

Notes regarding Nilsen’s MSE

Some side-notes that might be of use

Potentially relevant conferences and workshops

References consulted and annotated (Mendeley)

Update: 16 MAY 2018

GMSE v0.4.0.3 is now available on CRAN. A new website for GMSE has also been launched. This website was built with the R package pkgdown, recently released on CRAN. The site contains all of the vignettes and documentation for GMSE, and also includes a link to this lab notebook. A submission of the accompanying manuscript will soon be uploaded on bioRxiv.

Update: 14 MAY 2018

A new GMSE v0.4.0.3 has now been pushed to the master branch on GitHub and has been submitted to CRAN. The biggest update in this new version is a series of vignettes, plus a minor improvement to the genetic algorithm. More updates will follow soon, including some re-organisation of the GMSE project and a new manuscript submission.

Update: 13 APR 2018

I have re-worked the way that a manager restimates how the change in their policy affects users’ actions. The new new_act function in the genetic algorithm (games.c) performs well for getting more precise cost settings. The former way of doing it was much more of a blunt instrument, and it had a ceiling issue – that is, the manager would believe that higher costs caused fewer actions even when the resulting cost was over the users’ budgets.

/* =============================================================================
* This function updates an action based on the change in costs & paras
*     old_cost: The old cost of an action, as calculated in policy_to_counts
*     new_cost: The new cost of an action, as calculated in policy_to_counts
*     paras: Vector of global parameters
* ========================================================================== */
int new_act(double old_cost, double new_cost, double old_act, double *paras){

int total_acts;
double users, max_to_spend, acts_per_user, cost_per_user, total_cost;
double pr_on_act, budget_for_act, mgr_budget, min_cost;

users        = paras[54] - 1; /* Minus one for the manager */
min_cost     = paras[96];     /* Minimum cost of an action */
max_to_spend = paras[97];     /* Maximum per user budget   */
mgr_budget   = paras[105];    /* Manager's total budget    */

total_cost    = 0.0;
if(old_cost < mgr_budget){
total_cost    = old_act * old_cost; /* Total cost devoted to action */
}

cost_per_user = (total_cost / users); /* Cost devoted per user */
pr_on_act     = cost_per_user / max_to_spend;    /* Pr. devoted to action */

/* Assume that the proportion of the budget a user spends will not change */
budget_for_act = max_to_spend * pr_on_act;

/* Calculate how many actions to expect given acts per user and users */
acts_per_user  = budget_for_act / (new_cost + min_cost);
total_acts     = (double) users * acts_per_user;

return(total_acts);
}

This new way of assessing how users will act is now the function to be run in the background of all manager genetic agorithms. Very nicely, this also resolves an annoyance with the maximum allowed budgets. Previously, it was unclear why maximum budgets greater than 10000 were causing problems (managers were making bad predictions). I have now set the maximum budget to an order of magnitude higher, and there are no longer any apparent issues. A new version of GMSE will soon have this update.

Update: 21 DEC 2017

New Issue #40: Age distribution bump

Running simulations using gmse_apply, jeremycusack noticed a small but noticeable sharp decline in the population size at a generation equal to the maximum age of resources in the population (used a maximum age of 20). This decline is caused by the initial seed of resources having a uniform age distribution. In the first generation, these resources reproduce offspring that all have an age of zero, leading to an age structure in the population with many zero age individuals and a uniform distribution of ages greater than zero. The initial seed of individuals with random ages died gradually, but there were enough individuals in the initial offspring cohort that made it to the maximum age for it to have a noticeable effect in decreasing population size (i.e., all of these resources died on the maximum_age + 1 time step).

This effect can be avoided entirely given sufficient burn in generations of a model, and is less of a problem when the maximum age is low because this allows the age distribution to stabilise sooner. Further, using gmse_apply can avoid the issue by directly manipulating resources ages after the initial generation. Nevertheless, it would be useful to have a different default of age distributions instead of a uniform distribution.

One way to do this would be to find the age ($$A$$) at which a resource is expected to be alive with a probability of $$0.5$$, after accounting for mortality ($$\mu$$). This is simply calculated below:

$$(1 - \mu)^A = 0.5$$

The above can be re-arranged to find A,

$$A = \frac{log(0.5)}{log(1 - \mu)}$$.

Note that we could use a switch function (or something like it in R) to make $$A = 0$$ when $$\mu = 1$$, and revert to a uniform distribution of $$\mu = 0$$ (though this should rarely happen).

The value of $$\mu$$ would depend on res_death_type, and be processed in make_resource, which is used in both gmse and gmse_apply. If res_death_type = 1 (density independent, rarely used), then \mu is simply equal to remov_pr. If res_death_type = 2 (density dependent), then \mu could be found perhaps using something like the following:

mu = (RESOURCE_ini * lambda) / (RESOURCE_ini + RESOURCE_ini * lambda) gi This would get a value that is at least proportional to expected mortality rate of a resource (if res_death_type = 3, then we could use the some of types 1 and 2). Overall, the documentation should perhaps recommend finding a stable set of age distributions for a particular set of parameter combinations when using gmse_appy (i.e., through simulation), then using that distribution as an initial condition. But something like the above could probably get close to whatever the stable age distribution would be, at least close enough to make the decline in population size trivial.

I will start to consider some of the above as a potential default for the next version of GMSE. The best way to do this is probably to look at how code from the res_remove function in the resource.c file can best be integrated into a function called by the R function make_resource (i.e., either use the equations, or estimates of them, or somehow call res_remove directly).

Update: 13 DEC 2017

Improved convergence criteria

I have introduced and then immediately resolved Issue #39.

The convergence criteria has now been fixed with commit f598d8e52b47ef2017cac13d09aac1fb7aa6b506. To do this, I re-configured some of the genetic algorithm code into easier to read functions for checking the fitness increase. Now two separate ways of checking the increase in fitness from one genetic algorithm generation to the next exist; one for managers and one for users. This is needed because user fitness values are greater than zero and increase as their utility is maximised, but manager fitness values are less than zero and increase toward zero as their utility is maximised. The genetic algorithm now checks for a percentage improvement in fitness.

Now the default value of converge_crit equals 1, which means it does actually play a role sometimes (or is expected to). The genetic algorithm will continue until the percent increase in fitness from the previous generation is less than one percent. In practice, this doesn’t noticeably affect much, but it does allow better strategies to be found more quickly, and without having to play with ga_mingen to find them under extreme parameter settings (e.g., huge budgets and rapid shifts in abundance).

The new fix has now been checked and built with Winbuilder into v0.3.2.0, but I am leaving this on the development branch for now in anticipation of other potential improvements to be made soon.

Update: 1 NOV 2017

CRAN ready GMSE v0.3.1.7 – more flexibility, better error messages

I have now completed some substantial coding of error messages, which will be called in both gmse and gmse_apply. Essentially, these provide some help to software users who parameterise their models in a way that does not work with GMSE. For example, the parameter stakeholders is set to equal a negative number, an error message will be returned that informs the user that at least one stakeholder is required in the model. These error messages become a bit more important in gmse_apply, where it is possible for users to include arguments that don’t make sense (e.g., arrays of incorrect dimensions, or arguments that contradict one another).

The function gmse_apply has also been improved to make looping it easier. What had been happening during testing was that we were finding it all to easy to crash R by reading in parameters that contradicted one another (e.g., changing setting the landscape dimensions through land_dim_1 and land_dim_2 caused a crash when also trying to add in a LAND of different dimension – now this returns an error that LAND and land_dim_1 disagree about landscape size). This has been resolved in two ways. First, I have included many error messages meant to catch bad and contradictory arguments in gmse_apply (and, to a lesser extent gmse); it is still possible to crash R by setting things incorrectly, but you have to work very hard to do it – i.e., it almost has to be deliberate, as far as I can tell. Second, I have added the argument old_list to gmse_apply, which is FALSE by default, but can instead take the output of a previous full list return of gmse_apply (where get_res = Full). An element of the full list includes the basic output from which key parameters can be pulled. As a reminder, the basic gmse_apply output looks like the below.

$resource_results [1] 1062$observation_results
[1] 680.2721

$manager_results resource_type scaring culling castration feeding help_offspring policy_1 1 NA 110 NA NA NA$user_results
resource_type scaring culling castration feeding help_offspring tend_crops kill_crops
Manager             1      NA       0         NA      NA             NA         NA         NA
user_1              1      NA       9         NA      NA             NA         NA         NA
user_2              1      NA       9         NA      NA             NA         NA         NA
user_3              1      NA       9         NA      NA             NA         NA         NA
user_4              1      NA       9         NA      NA             NA         NA         NA

An example gmse_apply used in a loop is below.

to_scare <- FALSE;
sim_old  <- gmse_apply(scaring = to_scare, get_res = "Full", stakeholders = 6);
sim_sum  <- matrix(data = NA, nrow = 20, ncol = 7);
for(time_step in 1:20){
sim_new               <- gmse_apply(scaring = to_scare, get_res = "Full",
old_list = sim_old);
sim_sum[time_step, 1] <- time_step;
sim_sum[time_step, 2] <- sim_new$basic_output$resource_results[1];
sim_sum[time_step, 3] <- sim_new$basic_output$observation_results[1];
sim_sum[time_step, 4] <- sim_new$basic_output$manager_results[2];
sim_sum[time_step, 5] <- sim_new$basic_output$manager_results[3];
sim_sum[time_step, 6] <- sum(sim_new$basic_output$user_results[,2]);
sim_sum[time_step, 7] <- sum(sim_new$basic_output$user_results[,3]);
sim_old               <- sim_new;
print(time_step);
}
colnames(sim_sum) <- c("Time", "Pop_size", "Pop_est", "Scare_cost",
"Cull_cost", "Scare_count", "Cull_count");

The ouput sim_sum is shown below.

      Time Pop_size   Pop_est Scare_cost Cull_cost Scare_count Cull_count
[1,]    1      733  839.0023         NA       110          NA         54
[2,]    2      768  702.9478         NA       110          NA         54
[3,]    3      824  725.6236         NA       110          NA         54
[4,]    4      933  907.0295         NA       110          NA         54
[5,]    5     1180  816.3265         NA       110          NA         54
[6,]    6     1345 1224.4898         NA        10          NA        426
[7,]    7     1114 1269.8413         NA        10          NA        425
[8,]    8      820  884.3537         NA       110          NA         54
[9,]    9      952  793.6508         NA       110          NA         54
[10,]   10     1101  884.3537         NA       110          NA         54
[11,]   11     1299 1111.1111         NA        12          NA        402
[12,]   12     1079  907.0295         NA       110          NA         54
[13,]   13     1227 1564.6259         NA        10          NA        431
[14,]   14      934  839.0023         NA       110          NA         54
[15,]   15     1065 1133.7868         NA        10          NA        423
[16,]   16      768  725.6236         NA       110          NA         54
[17,]   17      869  929.7052         NA       110          NA         54
[18,]   18      949  907.0295         NA       110          NA         54
[19,]   19     1049  884.3537         NA       110          NA         54
[20,]   20     1200 1020.4082         NA        64          NA         90

We can take advantage of gmse_apply to dynamically change parameter values mid-loop. For example, below shows the same code, but with a policy of scaring introduced on time step 10.

to_scare <- FALSE;
sim_old  <- gmse_apply(scaring = to_scare, get_res = "Full", stakeholders = 6);
sim_sum  <- matrix(data = NA, nrow = 20, ncol = 7);
for(time_step in 1:20){
sim_new               <- gmse_apply(scaring = to_scare, get_res = "Full",
old_list = sim_old);
sim_sum[time_step, 1] <- time_step;
sim_sum[time_step, 2] <- sim_new$basic_output$resource_results[1];
sim_sum[time_step, 3] <- sim_new$basic_output$observation_results[1];
sim_sum[time_step, 4] <- sim_new$basic_output$manager_results[2];
sim_sum[time_step, 5] <- sim_new$basic_output$manager_results[3];
sim_sum[time_step, 6] <- sum(sim_new$basic_output$user_results[,2]);
sim_sum[time_step, 7] <- sum(sim_new$basic_output$user_results[,3]);
sim_old               <- sim_new;
if(time_step == 10){
to_scare <- TRUE;
}
print(time_step);
}
colnames(sim_sum) <- c("Time", "Pop_size", "Pop_est", "Scare_cost",
"Cull_cost", "Scare_count", "Cull_count");

The above simulation results in the following output for sim_sum.

      Time Pop_size   Pop_est Scare_cost Cull_cost Scare_count Cull_count
[1,]    1      745  657.5964         NA       110          NA         54
[2,]    2      805 1111.1111         NA        12          NA        400
[3,]    3      473  634.9206         NA       110          NA         54
[4,]    4      504  566.8934         NA       110          NA         54
[5,]    5      577  498.8662         NA       110          NA         54
[6,]    6      600  430.8390         NA       110          NA         54
[7,]    7      648  612.2449         NA       110          NA         54
[8,]    8      714  702.9478         NA       110          NA         54
[9,]    9      813  612.2449         NA       110          NA         54
[10,]   10      914 1020.4082         NA        64          NA         90
[11,]   11     1011 1179.1383         57        10          49        301
[12,]   12      858  725.6236         10       110         193         37
[13,]   13     1011 1043.0839         37        30           0        198
[14,]   14      989 1043.0839         57        30           0        198
[15,]   15      983 1065.7596         48        20          10        270
[16,]   16      851  839.0023         10       110         193         37
[17,]   17      962 1111.1111         38        12          58        306
[18,]   18      783  612.2449         10       110         193         37
[19,]   19      862  816.3265         10       110         193         37
[20,]   20      963  702.9478         10       110         182         38

Hence, in addition to all of the other benefits of gmse_apply, one new feature is that we can use it to study change in policy availability – in this case, what happens when scaring is introduced as a possible policy option. Similar things can be done, for example, to see how manager or user power changes over time. In the example below, users’ budgets increase by 100 every time step, with the manager’s budget remaining the same. The consequence appears to be decreased population stability and a higher likelihood of extinction.

ub          <- 500;
sim_old     <- gmse_apply(get_res = "Full", stakeholders = 6, user_budget = ub);
sim_sum     <- matrix(data = NA, nrow = 20, ncol = 6);
for(time_step in 1:20){
sim_new               <- gmse_apply(get_res = "Full", old_list = sim_old,
user_budget = ub);
sim_sum[time_step, 1] <- time_step;
sim_sum[time_step, 2] <- sim_new$basic_output$resource_results[1];
sim_sum[time_step, 3] <- sim_new$basic_output$observation_results[1];
sim_sum[time_step, 4] <- sim_new$basic_output$manager_results[3];
sim_sum[time_step, 5] <- sum(sim_new$basic_output$user_results[,3]);
sim_sum[time_step, 6] <- ub;
sim_old               <- sim_new;
ub                    <- ub + 100;
print(time_step);
}
colnames(sim_sum) <- c("Time", "Pop_size", "Pop_est", "Cull_cost", "Cull_count",
"User_budget");

The output of sim_sum is below.

      Time Pop_size   Pop_est Cull_cost Cull_count User_budget
[1,]    1     1215 1405.8957        10        292         500
[2,]    2     1065 1224.4898        10        336         600
[3,]    3      833  680.2721       110         36         700
[4,]    4      936  907.0295       110         42         800
[5,]    5     1174 1224.4898        10        401         900
[6,]    6      887  521.5420       110         54        1000
[7,]    7      988  680.2721       110         60        1100
[8,]    8     1084  975.0567       110         60        1200
[9,]    9     1208  861.6780       110         66        1300
[10,]   10     1360 1133.7868        10        520        1400
[11,]   11      975  861.6780       110         78        1500
[12,]   12     1079 1156.4626        10        560        1600
[13,]   13      597  770.9751       110         90        1700
[14,]   14      595  476.1905       110         96        1800
[15,]   15      586  612.2449       110        102        1900
[16,]   16      584  770.9751       110        108        2000
[17,]   17      557  589.5692       110        114        2100
[18,]   18      519  521.5420       110        120        2200
[19,]   19      469  521.5420       110        120        2300
[20,]   20      430  453.5147       110        126        2400

There is an important note to make about changing arguments to gmse_apply when old_list is being used: The function gmse_apply is trying to avoid a crash, so the function will accomodate parameter changes by rebuilding data structures if necessary. For example, if you change the number of stakeholders (and by including an argument stakeholders to gmse_apply, it is assumed that stakeholders are changing even they are not), then a new array of agents will need to be built. If you change landscape dimensions (or just include the argument land_dim_1 or land_dim_2), then a new landscape willl be built. This is mentioned in the documentation.

GMSE v0.3.3.7 passes all CRAN checks in Rstudio. I will make sure that the code works with win-builder, then prepare the new submission. Alternatively, as always, the newest GMSE version can be downloaded through GitHub if you have devtools installed in R.

devtools::install_github("bradduthie/GMSE")

I will soon update the manuscript for GMSE and upload it to biorXiv.

Update: 23 OCT 2017

Bug fix concerning density-based estimation

An error with density-based resource estimation (observe_type = 0) at very high values of agent_view was identified by Jeremy. When managers had a view of the landscape that encompassed a number of cells that was calculated to be larger than the actual number of landscape cells (as defined by land_dim_1 * land_dim_2), the manager would understimate actual population size. This occurred only in the manager.c file and not in the equivalent R function shown during plotting. The bug was fixed in commit a916b8f8a40041b5f08984cf73348108482dde59 with a simple if statement. This has therefore been resolved in a patched GMSE v0.3.1.3, which is now availabe on GitHub.

Update: 19 OCT 2017

Bug fix concerning resource movement

An error with the res_move_obs parameter was identified by Jeremy. This parameter was supposed to only affect resource movement during observation, but an if statement corrected in commit 5eeb88d285af57984171e7d72410659b3b441af3 was causing res_move_obs = FALSE to stop moving entirely in the resource model. This has now been resolved in a patched GMSE v0.3.1.1, which is now available on GitHub.

New option for removal of resources

A new option has been included for the argument res_death_type. By setting res_death_type = 3 in gmse or gmse_apply, resources can experience both density dependent (caused by res_death_K) and density independent (caused by remove_pr) removal simultaneously. Effects of each are independent of one another (i.e., both processes occur simultaneously, so the calculation of population size affecting removal due to carrying capacity includes resources that might experience density independent mortality).

Update: 16 OCT 2017

New group_think parameter in GMSE v0.3.1.0

A new group_think parameter has been developed by Jeremy and me, and included into an updated v0.3.1.0. This parameter is defined as FALSE be default, but when set to be TRUE will cause all users to act as a single block instead of independently. In the code, what happens is that a single user (user ID number 2) runs through the genetic algorithm, but then instead of having the resulting actions apply to this user, they apply to all users so that the genetic algorithm only needs to be run once in the user model. This decreases simulation time, particularly when there are a lot of users to model, but at a cost of removing all variation in actions among users. The group_think parameter can be defined in both gmse() and gmse_apply(), but I have not added it as an option in gmse_gui().

Update: 13 OCT 2017

GMSE v0.3.0.0 now available with gmse_apply

The gmse_apply function is now available on a new GMSE version 0.3.0.0. (minor tweaks to other functions have also been made, but nothing that changes the user experience of gmse – mostly typos corrected in the documentation). The new function allows software users to integrate their own submodels (resource, observation, manager, and user) into GMSE, or to use their own submodels entirely within a single function.

GMSE apply function

The gmse_apply function is a flexible function that allows for user-defined sub-functions calling resource, observation, manager, and user models. Where such models are not specified, GMSE submodels ‘resource’, ‘observation’, ‘manager’, and ‘user’ are run by default. Any type of sub-model (e.g., numerical, individual-based) is permitted as long as the input and output are appropriately specified. Only one time step is simulated per call to gmse_apply, so the function must be looped for simulation over time. Where model parameters are needed but not specified, defaults from gmse are used.

gmse_apply arguments

• res_mod The function specifying the resource model. By default, the individual-based resource model from gmse is called with default parameter values. User-defined functions must either return an unnamed matrix or vector, or return a named list in which one list element is named either ‘resource_array’ or ‘resource_vector’, and arrays must follow the format of GMSE in terms of column number and type (if there is only one resource type, then the model can also just return a scalar value).

• obs_mod The function specifying the observation model. By default, the individual-based observation model from gmse is called with default parameter values. User-defined functions must either return an unnamed matrix or vector, or return a named list in which one list element is named either ‘observation_array’ or ‘observation_vector’, and arrays must follow the format of GMSE in terms of column number and type (if there is only one resource type, then the model can also just return a scalar value).

• man_mod The function specifying the manager model. By default, the individual-based manager model that calls the genetic algorithm from gmse is used with default parameter values. User-defined functions must either return an unnamed matrix or vector, or return a named list in which one list element is named either ‘manager_array’ or ‘manager_vector’, and arrays must follow the (3 dimensional) format of the ‘COST’ array in GMSE in terms of column numbers and types, with appropriate rows for interactions and layers for agents (see documentation of GMSE for constructing these, if desired). User defined manager outputs will be recognised as costs by the default user model in gmse, but can be interpreted differently (e.g., total allowable catch) if specifying a custom user model.

• use_mod The function specifying the user model. By default, the individual-based user model that calls the genetic algorithm from gmse is used with default parameter values. User-defined functions must either return an unnamed matrix or vector, or return a named list in which one list element is named either ‘user_array’ or ‘user_vector’, and arrays must follow the (3 dimensional) format of the ‘ACTION’ array in GMSE in terms of column numbers and types, with appropriate rows for interactions and layers for agents (see documentation of GMSE for constructing these, if desired).

• get_res How the output should be organised. The default ‘basic’ attempts to distill results down to their key values from submodel outputs, including resource abundances and estimates, and manager policy and actions. An option ‘custom’ simply returns a large list that includes the output of every submodel. Any other option (e.g. ‘full’) will return a massive list with all of the input, output, and parameters used to run gmse_apply.

• Arguments passed to user-defined functions, and passed to modify default parameter values that would otherwise be called for gmse default models. Any argument that can be passed to gmse can be specified explicitly, just as if it were an argument to gmse. Similarly, any argument taken by a user-defined function should be specified, though the function will work if the user-defined function has a default that is not specified explicitly.

Example uses of gmse_apply

A simple run of gmse_apply() will return one generation of gmse using default submodels and parameter values.

sim <- gmse_apply();

For sim, the default ‘basic’ results are returned as below.

$resource_results [1] 1102$observation_results
[1] 1179.138

$manager_results scaring culling castration feeding help_offspring policy NA 10 NA NA NA$user_results
resource_type scaring culling castration feeding help_offspring tend_crops kill_crops
Manager             1      NA       0         NA      NA             NA         NA         NA
user_2              1      NA      70         NA      NA             NA         NA         NA
user_3              1      NA      75         NA      NA             NA         NA         NA
user_4              1      NA      69         NA      NA             NA         NA         NA
user_5              1      NA      74         NA      NA             NA         NA         NA

Note in the case above we have the total abundance of resources returned, the estimate of resource abundance from the observation function, the costs the manager sets for the only available action of culling, and the number of culls attempted by each user.

The above was produced by all of the individual-based functions that are default in GMSE; custom generated subfunctions can instead be included provided that they fit the specifications described above. For example, we can define a very simple logistic growth function to send to res_mod instead.

alt_res <- function(X, K = 2000, rate = 1){
X_1 <- X + rate*X*(1 - X/K);
return(X_1);
}

The above function takes in a population size of X and returns a value X_1 based on the population intrinsic growth rate rate and carrying capacity K. Iterating the logistic growth model by itself under default parameter values with a starting population of 100 will cause the population to increase to carrying capacity in roughly 7 generations. The function can be substituted into gmse_apply to use it instead of the default GMSE resource model.

sim <- gmse_apply(res_mod = alt_res, X = 100, rate = 0.3);

The gmse_apply function will find the parameters it needs to run the alt_res function in place of the default resource function, either by running the default function values (e.g., K = 2000) or values specified directly into gmse_apply (e.g., X = 100 and rate = 0.3). If an argument to a custom function is required but not provided either as a default or specified in gmse_apply, then an error will be returned.

To integrate across different types of submodels, gmse_apply translates between vectors and arrays between each submodel. For example, because the default GMSE observation model requires a resource array with particular requirements for column identites, when a resource model subfunction returns a vector, or a list with a named element ‘resource_vector’, this vector is translated into an array that can be used by the observation model. Specifically, each element of the vector identifies the abundance of a resource type (and hence will usually be just a single value denoting abundance of the only focal population). If this is all the information provided, then a resource_array will be made with default GMSE parameter values with an identical number of rows to the abundance value (floored if the value is a non-integer; non-default values can also be put into this transformation from vector to array if they are specified in gmse_apply, e.g., through an argument such as lambda = 0.8). Similarly, a resource_array is also translated into a vector after the default individual-based resource model is run, should the observation model require simple abundances instead of an array. The same is true of observation_vector and observation_array objects returned by observation models, of manager_vector and manager_array (i.e., COST) objects returned by manager models, and of user_vector and user_array (i.e., ACTION) objects returned by user models. At each step, a translation between the two is made, with necessary adjustments that can be tweaked through arguments to gmse_apply when needed. Alternative observation, manager, and user, submodels, for example, are defined below; note that each requires a vector from the preceding model.

# Alternative observation submodel
alt_obs <- function(resource_vector){
X_obs <- resource_vector - 0.1 * resource_vector;
return(X_obs);
}

# Alternative manager submodel
alt_man <- function(observation_vector){
policy <- observation_vector - 1000;
if(policy < 0){
policy <- 0;
}
return(policy);
}

# Alternative user submodel
alt_usr <- function(manager_vector){
harvest <- manager_vector + manager_vector * 0.1;
return(harvest);
}

All of these submodels are completely deterministic, so when run with the same parameter combinations, they produce replicable outputs.

gmse_apply(res_mod = alt_res, obs_mod = alt_obs,
man_mod = alt_man, use_mod = alt_usr, X = 1000);

The above, for example, produces the following output (Note that the X argument needs to be specified, but the rest of the subfunctions take vectors that gmse_apply recognises will become available after a previous submodel is run).

$resource_results [1] 1500$observation_results
[1] 1350

$manager_results [1] 350$user_results
[1] 385

Note that the manager_results and user_results are ambiguous here, and can be interpreted as desired – e.g., as total allowable catch and catches made, or as something like costs of catching set by the manager and effort to catching made by the user. Hence while manger output is set in terms of costs of performing each action, and user output is set in terms of action attempts, this need not be the case when using gmse_apply (though it should be recognised when using default GMSE manager and user functions).

GMSE default submodels can be added in at any point.

gmse_apply(res_mod = alt_res, obs_mod = observation,
man_mod = alt_man, use_mod = alt_usr, X = 1000)

The above produces the results below.

$resource_results [1] 1500$observation_results
[1] 1655.329

$manager_results [1] 655.3288$user_results
[1] 720.8617

If we wanted to, for example, specify a simple resource and observation model, but then take advantage of the genetic algorithm to predict policy decisions and user actions, we could use the default GMSE manager and user functions (written below explicitly, though this is not necessary).

gmse_apply(res_mod = alt_res, obs_mod = alt_obs,
man_mod = manager, use_mod = user, X = 1000)

The above produces the output below returning culling costs and culling actions attempted by four users (note that the default manager target abundance is 1000).

$resource_results [1] 1500$observation_results
[1] 1350

$manager_results scaring culling castration feeding help_offspring policy NA 10 NA NA NA$user_results
resource_type scaring culling castration feeding help_offspring tend_crops kill_crops
Manager             1      NA       0         NA      NA             NA         NA         NA
user_2              1      NA      70         NA      NA             NA         NA         NA
user_3              1      NA      70         NA      NA             NA         NA         NA
user_4              1      NA      71         NA      NA             NA         NA         NA
user_5              1      NA      73         NA      NA             NA         NA         NA

Instead of using the gmse function, we might simulate multiple generations by calling gmse_apply through a loop, reassigning outputs where necessary for the next generation (where outputs are not reassigned, new defaults will be inserted in their place, so, e.g., if we were to just loop without reassigning any variables, nothing would update and we would be running the same model, effectively, multiple times). Below shows how this might be done.

sim1      <- gmse_apply(get_res = "full", lambda = 0.3);
RESOURCES <- sim1$resource_array; LAND <- sim1$LAND;
PARAS     <- sim1$PARAS; results <- matrix(dat = NA, nrow = 40, ncol = 4); for(time_step in 1:40){ sim_new <- gmse_apply(RESOURCES = RESOURCES, LAND = LAND, PARAS = PARAS, COST = COST, ACTION = ACTION, stakeholders = 10, get_res = "full", agent_view = 20); results[time_step, 1] <- sim_new$resource_vector;
results[time_step, 2] <- sim_new$observation_vector; results[time_step, 3] <- sim_new$manager_vector;
results[time_step, 4] <- sim_new$user_vector; RESOURCES <- sim_new$resource_array;
LAND      <- sim_new$LAND; PARAS <- sim_new$PARAS;
COST      <- sim_new$COST; ACTION <- sim_new$ACTION;
}

colnames(results) <- c("Abundance", "Estimate", "Cull_cost", "Cull_attempts");

The above results in the following output for results.

      Abundance  Estimate Cull_cost Cull_attempts
[1,]      1195 1165.9726        10   716
[2,]      1045  939.9167       110   461
[3,]      1160 1160.0238        10   715
[4,]      1056 1183.8192        10   715
[5,]      1014  850.6841       110   468
[6,]      1171 1237.3587        10   717
[7,]      1026  993.4563       110   464
[8,]      1202  957.7632       110   464
[9,]      1394 1469.3635        10   702
[10,]      1333 1457.4658        10   702
[11,]      1277 1397.9774        10   702
[12,]      1175 1415.8239        10   702
[13,]      1088  701.9631       110   468
[14,]      1275 1207.6145        10   718
[15,]      1200 1332.5402        10   718
[16,]      1116 1029.1493        45   512
[17,]      1249 1814.3962        10   699
[18,]      1141 1273.0518        10   722
[19,]      1019  963.7121       110   455
[20,]      1216 1629.9822        10   708
[21,]      1088 1130.2796        10   708
[22,]       988 1035.0982        38   537
[23,]      1056 1029.1493        45   505
[24,]      1154  749.5538       110   463
[25,]      1344 1499.1077        10   722
[26,]      1268 1386.0797        10   712
[27,]      1165 1493.1588        10   707
[28,]      1061 1070.7912        19   633
[29,]      1019 1076.7400        17   663
[30,]       961  600.8328       110   457
[31,]      1135  874.4795       110   450
[32,]      1338 1189.7680        10   701
[33,]      1275 1600.2380        10   710
[34,]      1174 1362.2844        10   709
[35,]      1104 1112.4331        12   685
[36,]      1003 1302.7960        10   715
[37,]       828 1183.8192        10   712
[38,]       649  785.2469       110   462
[39,]       739 1023.2005        56   488
[40,]       813  910.1725       110   455

Note that managers increase the cost of culling based on the time step’s estimated abundance, and user culling attempts decrease when culling costs increase.

In addition to the flexibility of allowing user-defined submodels, gmse_apply is also useful for modellers who might be interested in simulating processes not currently available in gmse by itself. For example, if we wanted to model a sudden environmental perturbation decreasing population size, or a sudden influx of new users, after 30 generations, we could do so in the loop.

In the near future, the gmse_apply function will be included in the GMSE vignette and submitted to CRAN with the rest of v0.3.0.0 – in the mean time, I believe that all major bugs have been ironed out, but please let me know or report an issue if you are able to crash the function (i.e., if you run it and it causes R to crash – you should always get an error message before this happens).

To download the latest GMSE v0.3.0.0, simply run the below in R (make sure that devtools is installed).

devtools::install_github("bradduthie/GMSE")

I welcome any feedback, and I expect to submit an update to CRAN around late October.

Update: 12 OCT 2017

New function gmse_apply complete and tested

I have now completed the gmse_apply function, which exploits the full modularity of GMSE by allowing software users to develop their own sub-functions and string them together with any combination of GMSE default sub-functions. As a brief summary, gmse_apply includes the following features:

• Any coherent function can be used for a resource model as long as it returns an unlabelled vector or matrix, or returns a list with at least one element labelled ‘resource_vector’ or ‘resource_array’.
• Similarly, any function can be used for a resource model as long as it takes an argument called ‘resource_vector’ or ‘resource_array’ and returns an unnamed vector or matrix, or a list with at least one element labelled ‘observation_vector’ or ‘observation_array’.
• Any function works for the manager model if it takes an argument called ‘observation_vector’ or ‘observation_array’ and returns an unlabelled vector or matrix, or a list with at least one element labelled ‘manager_vector’ or ‘manager_array’.
• Any function works for the user model if it takes an argument called ‘manager_vector’ or ‘manager_array’ and returns an unlabelled vector or matrix, or a list with at least one element labelled ‘user_vector’ or ‘user_array’.
• Returns one of three options, including (1) a breakdown of the key results (default), (2) all of the model output, or (3) a list of just the sub-model outputs.

Any arguments for custom user functions can simply be passed along by specifying them in gmse_apply. For example, if we have a custom resource function alt_res below:

alt_res <- function(X = 1000, K = 2000, r = 1){
X_1 <- X + r*X*(1 - X/K);
return(X_1);
}

We can simply include the above in gmse_apply as follows to use the very simple logistic growth sub-model with the individual-based submodels that are defaults of GMSE.

sim_app <- gmse_apply(res_mod = alt_res);

The gmse_apply function simply adds in GMSE defaults for unspecified models, but we can specify them too.

sim_app <- gmse_apply(res_mod = alt_res, obs_mod = observation);

To adjust parameters in the alternative resource model, simply add in the arguments as below.

sim_app <- gmse_apply(res_mod = alt_res, X = 2000, K = 5000, r = 1.2);

The gmse_apply function will know where to place them, and update them should they be needed for other models.

I will give a more lengthy description of how to use gmse_apply tomorrow, when I push GMSE v0.3.0.0 to the master branch of GitHub and advertise the update.

Update: 6 OCT 2017

Compensation suggestion

A suggestion from Jeremy to include a compensation option for users. Users could devote some of their budget to compensation, then managers could compensate a proportion of their damaged yield. Implementing this will require consideration from the manager’s perspective with respect to the genetic algorithm – the users’ perspective will be easier because a user can remember their previous losses and assess compensation versus culling. Managers might have to think about how compensation could incentivise non-culling, but this might actually already work given the way the manager anticipates actions; more investigation into this will be useful following the finalisation of gmse_apply(), which is in progress.

Update: 28 SEP 2017

Progress has been made on the gmse_apply() function. My goal is to make this as modular as possible – to allow any four functions to be included in the GMSE framework, including arbitrary arguments to each function. The gmse_apply() function will recognise which arguments go along with which functions, and naturally string together results from one sub-function to the input of the next sub-function (though this will demand that the output from functions is labelled in a way that matches the arguments of the next function; e.g., if you have a ‘N_total’ as input for the observation model, then ‘N_total’ will either need to be labelled output of the resource model or specificied directly in gmse_apply()). Default submodels will be the IBMs used in gmse(), and where arguments are not specified by the software user in gmse_apply() (e.g., LAND) they will be built from default gmse() parameters.

Update: 27 SEP 2017

The GMSE GUI has been updated with all of the new features in version 0.2.3.0. The gmse_gui() function is likewise updated in a new patch version 0.2.3.1. I did this quickly because the GUI was actually easy to update; plans for the gmse_apply function are now also clear, and I hope to have a working function and version 0.3.0.0 by the end of the week, or by early next week.

Update: 26 SEP 2017

GMSE Version 0.2.3.0 on GitHub

I have pushed a new version 0.2.3.0 of GMSE onto the master branch of GitHub, which means that the most up-to-date version can be installed using the code below (make sure the devtools library is installed).

devtools::install_github("bradduthie/GMSE")

The new version includes multiple new features:

• A bug fix (see below) that was causing confusion between culling and castration actions.
• A new plot_gmse_effort function, which shows the conflict between manager targets and user actions more directly (see plots from 23 AUG 2017 notes).
• A new gmse_summary function, which takes the large output produced from gmse and returns a much easier to understand set of tables.
• A revised gmse_gui function that has better defaults and parameter organisation, which has been uploaded to shiny for use in a browser.
• Slightly adjusted default parameter values.
• Documentation for new functions, and error fixes in the old documentation.

To run a simple default simulation, the gmse function remains unchanged.

sim <- gmse();

To plot the effort of managers and users, use the below.

plot_gmse_effort(agents = sim$agents, paras = sim$paras,
ACTION = sim$action, COST = sim$cost);

Below summarises the results more cleanly, extracting key information from sim.

gmse_summary(sim);

And as before, the GUI can be called directly from the R console.

gmse_gui();

The GUI does not yet allow you to get a vew of the plot_gmse_effort output, or a gmse_summary, but this will be a goal for future versions of GMSE.

If able, I recommend updating to version 0.2.3.0 as soon as possible. In the coming few days, I will also add the gmse_apply function, primarily for developers who will benefit from a more modular way of using GMSE, allowing for different types of submodules to be used within the broader GMSE framework. When the new apply function has been added (and possibly the GUI improved), I will submit a new version 0.3.x.x to CRAN.

Bug Fix and tweaks to agent prediction

I have now fixed a bug in the code that was causing confusion between culling and castration. After recompiling and running simulations, manager and user actions improve. I have also made some minor changes to default gmse() options. Regarding the predicted consequences of manager and user actions (i.e., the predictions from the agents’ perspective that guid their decision making), I have adjusted some things to make them more in line with what is expected in the simulation as follows (recall that managers are interested in global abundance and users are interested specifically in how abundance affects themselves):

1. Scaring: Managers predict no change in resource abundance, while users predict a decrease of 1
2. Culling: Managers and users predict a decrease of $$1 + \lambda$$ (Note that this brings in knowledge of birth rate a priori – might want to allow for a change in this in the simulation, but it also seems realistic for agents to recognise that adults can reproduce, and a value is needed to reflect this)
3. Castration: Managers and users predict a decrease of $$\lambda$$
4. Feeding: Managers and users predict an increase of $$\lambda$$
5. Help offspring: Managers and users predict an increase of 1

These values are a bit more in line with what will actually happen, so we assume that managers and users are a bit more informed now. It also allows for a bit more differentiation among actions. Overall, the model appears to perform better now – meaning that managers and users appear to be better predictors of the conseuqneces of their actions.

Before finishing the gmse_apply() function, I will push an updated version of GMSE to GitHub with these changes, plus new plotting options.

Update: 25 SEP 2017

I have written a gmse_summary function (see below), which returns a simplified list that includes four elements, each of which is a table of data: 1. resources, a table showing time step in the first column, followed by resource abundance in the second column. 2. observations, a table showing time step in the first column, followed by the estimate of population size (produced by the manager) in the second column. 3. costs, a table showing time step in the first column, manager number in the second column (should always be zero), followed by the costs of each action set by the manager (policy); the far-right column indicates budget that is unused and therefore not allocated to any policy. 4. actions, a table showing time step in the first column, user number in the second column, followed by the actions of each user in the time step; additional columns indicate unused actions, crop yield on the user’s land (if applicable), and the number of resources that a user successfully harvests (i.e., ‘culls’).

At the moment, I have not added in the actual number of resources that a user culls. This will be added shortly, after which I will post a new function. Doing so is a bit more complicated because it requires me to go into the C code and make a recording every time it happens (see how I plan to do this below the function).

gmse_summary <- function(gmse_results){
time_steps <- dim(gmse_results$paras)[1]; parameters <- gmse_results$paras[1,];
#--- First get the resource abundances
res_types    <- unique(gmse_results$resource[[1]][,2]); resources <- matrix(dat = 0, nrow = time_steps, ncol = length(res_types) + 1); res_colna <- rep(x = NA, times = dim(resources)[2]); res_colna[1] <- "time_step"; for(i in 1:length(res_types)){ res_colna[i+1] <- paste("type_", res_types[i], sep = ""); } colnames(resources) <- res_colna; #--- Next get estimates abd the costs set by the manager observations <- matrix(dat = 0, nrow = time_steps, ncol = length(res_types) + 1); costs <- matrix(dat = NA, nrow = time_steps*length(res_types), ncol = 10); agents <- gmse_results$agents[[1]];
users   <- agents[agents[,2] > 0, 1];
actions <- matrix(dat  = NA, ncol = 13,
nrow = time_steps * length(res_types) * length(users));
c_row  <- 1;
a_row  <- 1;
for(i in 1:time_steps){
the_res            <- gmse_results$resource[[i]][,2]; manager_acts <- gmse_results$action[[i]][,,1];
resources[i, 1]    <- i;
observations[i, 1] <- i;
land_prod          <- gmse_results$land[[i]][,,2]; land_own <- gmse_results$land[[i]][,,3];
for(j in 1:length(res_types)){
#---- Resource abundance below
resources[i,j+1] <- sum(the_res == res_types[j]);
#---- Manager estimates below
target_row <- which(manager_acts[,1] == -2 &
manager_acts[,2] == res_types[j]);
estim_row  <- which(manager_acts[,1] ==  1 &
manager_acts[,2] == res_types[j]);
target <- manager_acts[target_row, 5];
#---- Cost setting below
costs[c_row, 1]  <- i;
costs[c_row, 2]  <- res_types[j];
estim_row    <- which(manager_acts[,1] ==  1 &
manager_acts[,2] == res_types[j]);
if(parameters[89] == TRUE){
costs[c_row, 3] <- manager_acts[estim_row,  8];
}
if(parameters[90] == TRUE){
costs[c_row, 4] <- manager_acts[estim_row,  9];
}
if(parameters[91] == TRUE){
costs[c_row, 5] <- manager_acts[estim_row,  10];
}
if(parameters[92] == TRUE){
costs[c_row, 6] <- manager_acts[estim_row,  11];
}
if(parameters[93] == TRUE){
costs[c_row, 7] <- manager_acts[estim_row,  12];
}
if(parameters[94] == TRUE){
costs[c_row, 8] <- parameters[97];
}
if(parameters[95] == TRUE){
costs[c_row, 9] <- parameters[97];
}
costs[c_row, 10] <- manager_acts[estim_row, 13] - parameters[97];
c_row <- c_row + 1;
#--- Action setting below
for(k in 1:length(users)){
usr_acts <- gmse_results$action[[i]][,,users[k]]; actions[a_row, 1] <- i; actions[a_row, 2] <- users[k]; actions[a_row, 3] <- res_types[j]; res_row <- which(usr_acts[,1] == -2 & usr_acts[,2] == res_types[j]); if(parameters[89] == TRUE){ actions[a_row, 4] <- usr_acts[res_row, 8]; } if(parameters[90] == TRUE){ actions[a_row, 5] <- usr_acts[res_row, 9]; } if(parameters[91] == TRUE){ actions[a_row, 6] <- usr_acts[res_row, 10]; } if(parameters[92] == TRUE){ actions[a_row, 7] <- usr_acts[res_row, 11]; } if(parameters[93] == TRUE){ actions[a_row, 8] <- usr_acts[res_row, 12]; } if(j == length(res_types)){ if(parameters[104] > 0){ land_row <- which(usr_acts[,1] == -1); if(parameters[95] > 0){ actions[a_row, 9] <- usr_acts[land_row, 10]; } if(parameters[94] > 0){ actions[a_row, 10] <- usr_acts[land_row, 11]; } } actions[a_row, 11] <- sum(usr_acts[, 13]); } if(parameters[104] > 0){ max_yield <- sum(land_own == users[k]); usr_yield <- sum(land_prod[land_own == users[k]]); actions[a_row, 12] <- 100 * (usr_yield / max_yield); } a_row <- a_row + 1; } } } cost_col <- c("time_step", "resource_type", "scaring", "culling", "castration", "feeding", "helping", "tend_crop", "kill_crop", "unused"); colnames(costs) <- cost_col; colnames(resources) <- res_colna; colnames(observations) <- res_colna; action_col <- c("time_step", "user_ID", "resource_type", "scaring", "culling", "castration", "feeding", "helping", "tend_crop", "kill_crop", "unused", "crop_yield", "harvested"); colnames(actions) <- action_col; the_summary <- list(resources = resources, observations = observations, costs = costs, actions = actions); return(the_summary); } To record kills, I think that the best way is to use the resource mortality adjustment column (at the moment, column 17 in C and 18 in R of the resource array). Mortality as of now is just adjusted to 1 in the event of a kill, and mortality occurs whenever a random probability is greater than or equal to 1. Hence, I can replace the 1 value with the user’s ID (for non-managers, this must be at least 1), and then the resource array will record the ID of the user that killed it at the particular time step. Note that this cannot be done for other adjustments such as growth rate or offspring production because the values are not interpreted as probabilities. I will do the above tomorrow, which should not take too long. I will then continue work on the gmse_apply function. Update: 22 SEP 2017 Currently, the gmse() function returns a list that includes all of the data produced by the model, some details of which are required for plotting. sim_results <- list(resource = RESOURCE_REC, observation = OBSERVATION_REC, paras = PARAS_REC, land = LANDSCAPE_REC, time_taken = total_time, agents = AGENT_REC, cost = COST_REC, action = ACTION_REC ); I think that this list is fine, perhaps necessary, to keep, but the ConFooBio group has also concluded that there should be some easier to understand summary of the data. I propose that some function written, a gmse_summary(), that summarises the results in an easier to understand way would be useful. The function could just be run as below. sim <- gmse(); sim_summary <- gmse_summary(sim); The output of gmse_summary() should be a list of all of the relevant information that a user might want to plot or analyse. It should include the following list elements. 1. sim_summary$resources
2. sim_summary$observations 3. sim_summary$costs
4. sim_summary\$actions

More might be needed, but the above should be a good starting point that will provide four clear data tables for the user. The tables will look like the below.

1. Resource abundances over time

time_step abundance
1 100
2 104
99 116
100 108

In the above, only the resource abundance is reported to the software user, though it might also be useful to have additional columns as well eventually.

2. Observation estimates of abundance over time

time_step estimated_abundance
1 102
2 101
99 121
100 112

In the above, only the estimate from the observaiton submodel is reported to the software user. Additional columns might also be useful for things like confidence intervals, though for now I’m not sure if this is needed.

3. Costs set in each time step

time_step manager scaring castration culling feeding helping unused
1 0 40 NA 60 NA NA 0
2 0 36 NA 62 NA NA 2
99 0 0 NA 100 NA NA 0
100 0 3 NA 97 NA NA 0

In the above, the manager number is always 0 because this is the number of the agent that has that role in GMSE. All impossible actions (specificed by the simulation) are labelled NA, while the possible scaring and culling actions are given values that correspond to the cost of each action for users in each time step. Hence the table summarises policy for each time step in a way that software users can interpret more cleanly.

4. Actions in each time step

time_step user scaring castration culling feeding helping tend_crop kill_crop unused crop_yield harvested
1 1 50 NA 50 NA NA NA NA 0 90 12
1 2 59 NA 40 NA NA NA NA 1 92 9
1 3 100 NA 0 NA NA NA NA 0 89 0
2 1 44 NA 66 NA NA NA NA 0 88 16
2 2 52 NA 48 NA NA NA NA 0 94 12
2 3 98 NA 0 NA NA NA NA 2 90 0
99 1 36 NA 63 NA NA NA NA 1 79 20
99 2 40 NA 60 NA NA NA NA 0 83 18
99 3 28 NA 72 NA NA NA NA 0 88 12
100 1 35 NA 62 NA NA NA NA 3 82 18
100 2 37 NA 63 NA NA NA NA 0 84 22
100 3 23 NA 77 NA NA NA NA 0 84 13

The above action table has more rows in it than the cost table because a row is needed for each user in each time step. This gives the software user full access to each individual user’s actions, and their results. Note that as above, castration, feeding, and helping, are not options. Additionally, in this hypothetical simulation, tending or killing crops are not options, so no actions are performed. Users divide their budget between scaring and culling in each time step. The last two columns also give useful information to the software user. The first is crop yield on the user’s owned land (should probably be NA if land_ownership = FALSE), which will reflect the percentage of the total possible yield (or maybe raw yield?) for each user – hencing allowing the table to direclty correlate actions with yield. The last column is the number of resources ‘harvested’, which I think should count successful ‘kills’ (rather than just actions devoted to culling). The realised culling might be lower than the actions devoted to culling, for example, if not enough resources are actually on the user’s land to cull. Additional statistics for each user could be added in as columns, but this seems a good place to start. This gmse_summary producing a list of the above four tables will be included in the next version of GMSE, along with the new plotting function highlighting the conflict itself, and the gmse_apply function discussed on 6 SEPT.

Update: 15 SEP 2017

Continued progress has been made on slides for an upcoming talk.

Update: 13 SEP 2017

I will be giving a talk on 19 September 2017 for the Mathematics and Statistics Group at the University of Stirling on GMSE as a general tool for management strategy evaluation. Slides for this talk will be available on GitHub.

Update: 8 SEP 2017

The alternative approach from Wednesday is being implemented smoothly. Passing user-defined functions in a modular way is possible, but inputs and outputs need to be carefully considered within gmse_apply(). The objective is to make things as easy and flexible as possible for the user, while also making sure that the function runs efficiently.

Update: 6 SEP 2017

A modular function for modellers

I am beginning work on a gmse_apply() function, which will improve the modularity of GMSE for developers. The goal behind this function will be to provide a simple tool for allowing developers to use their own resource and observation models and, with the correct inputs, take advantage of the manager and user functions. Hence, simple resource and observation models will be possible, but the flexibility of GMSE should be retained as much as possible. A few starting points include the following:

• The gmse_apply() function will run a single cycle of GMSE instead of multiple time steps.
• The function will take user-supplied functions as input for each sub-model.
• All relevant gmse() options will be acceptable, but not obvious, appearing as ... passes that the user can add if they want to change things. Otherwise, defaults will be used.
• The function will allow for both numerical and individual-based models, but numerical models will be the default.
• Numerical models will work with a vector instead of a table of resources. For now, this vector will include just a single estimate that is the actual resource population size (could be expanded to different resource types). A second vector will include estimates, produced by the observation model. The default manager model will take these vectors and run the genetic algorithm accordingly, as will the user model.
• There should be room for software users to supply their own manager and user models, but with the appropriate checks to ensure model compatibility (an early switch might organise things clearly).

Inputs and outputs of different functions will then include the following:

• Resource model
• Numerical will require a vector of population size values from the larger function (for now, this will almost always just be a single element of the population size). Harvesting will be output at the end, so that users can just add in a simple model (e.g., logistic growth) and not worry about adding the harvesting in here (as it does in MSEtools). LIkewise, the output will only be a vector of population size values, so very simple models will be possible. Additional model-specific parameters will be specified in the function itself. In other words, for a user inputting their own logistic growth model into gmse_apply(), the only thing required of the user-defined function will be the population size vector, and other parameters will be specified in the logistic function, e.g.: gmse_apply(resource = LV(pop, K = 100, r = 1, ...), observation = ..., type = "Numerical").
• Individual-based will default to the resource() function, though some options to input landscape and starting conditions will be needed (though these should also switch to default).
• Observation model
• Numerical will require a vector of population size values from the larger function, and a vector will be returned that is the observed estimate of the population size. Other options might affect the observation function itself, but no other input or output should be necessary, as in the resource function.
• Individual-based will default to the observation() function, as as with the resource model.
• Manager model
• Numerical will need to be tweaked so that it can accept direct estimates from the observation model, so the estimate_abundances function in manager.c can be bipassed entirely (i.e., I still think that we’ll want to call the c function as normal, but add an option). This can be arranged simply by reading in the observation vector (might have to re-structure to an array a bit better?) as the OBSERVATION array in c, but then instead of running estimate_abundances, allow an if-else statement to read this array as abun_est given a new value in paras, after which the genetic algorithm and set_action_costs can be run as normal. The irrelevant output can just be ignored by the user model in gmse_apply.
• Individual-based will default to the manager() function in its normal form. As with the others, this will require some eventual decisions about initialisation, but I can worry about this later.
• There will be an option built into gmse_apply() that would allow users to specify their own manager functions, but for the moment this is just going to be hidden because the genetic algorithm requires the COST and ACTION arrays to be in the correct form for use. Hence, if someone later wants to apply their own manager or user function, they will either need to get the input-output correct (at the moment) or (eventually) use some different input-output data structure that I make up later; but at some point, things would just collapse down to MSEtools. Or, rather, gmse_apply() just becomes a trivial function that includes four lines of compatible code, calling each model.
• User model
• Numerical will also need to be tweaked to return actions from the user model, but this shouldn’t be too difficult. A key here is that given these very simple models, the only real policy option is ‘culling’. We really don’t even need to call the user function if we just assume that users always cull as much as they possibly can. Nevertheless, I’ll call it and the genetic algorithm along with it, sticking in place the relevant parameters for a spatially implicit model. As with the manager model, there will need to be some sort of if-else where landscape based actions do not apply, and actions will be not applied to the resource array. Instead, the ACTION output from user() will just be translated in R to adjusting the vector.
• Individual-based will default to the user() function in its normal form. As with the others, this will require some eventual decisions about initialisation to be determined later.

As an alternative, at least to the implementation, I think that the call could be made at the level of the individual resource() and observation() R functions. This was kind of always the plan, but there’s a semi-dirty way to mix numerical resource and observation models with the full individual-based manager and user models. This can be done by adding a model option to be user defined through an if( is.function(model) == TRUE ) in the resource.R function. If the condition is satisfied, then resource() will shift to the user generated model. This can actually be done for all of the submodels very easily.

This alternative might be a better way to go. The aforementioned ‘dirty’ part of the technique might be to check to see if the output is in the correct form, then, if only a vector is returned – turn it into the correct form by making a data frame that has the same number of rows by calling make_resource. The type 1 values could correspond to vector elements. Admittedly, this could get slow for huge population sizes, but population sizes would have to be massive for R to slow down from simplying making a matrix with a lot of rows. In any case, it would at least standardise the input and output for the user of gmse_apply in a way that plays nice with everything else in GMSE.

Similarly, the observation function could also call make_resource if a vector is returned (since individual variation wouldn’t be relevant in the numerical model).

With this alternative approach, no changes to the C code need to be made – the inputs and outputs just need to be tweaked into a standardised way when a vector or scalar is returned from any user-defined model (small detail – population size needs to be an integer). This can be an option later for the user and manager models – though I’m not sure how this would work, exactly. A benefit here is that some parts of the model could concievably individual-based, with others being numerical – the trade-off being the requirement for discrete resource numbers and a very small amount of slowdown (which will almost certainly not be noticable for any resaonable model).

The gmse_apply function would then initialise a very small number of agents, and a small landscape (unless otherwise specified) in every run. The possibility of passing more options could be applied with a simple .... This would also require a sub-function build_para_vec, which would be used for the sole purpose of taking the list of options included (same as in gmse()) and passing it to the sub-function, with any functions not passed being assumed defaults (and most would be irrelevant). So the default function should then look like

gmse_apply(resource_function = "IBM", observation_function = "IBM", manager_function = "IBM", user_function = "IBM", res_number = 100, ...);

I think at least an initial population needs to be specified, but everything else can be left up to the user, with the elipses passing to the function building the parameter vector (which can also be called by gmse(), replacing some clutter). Overall, the function will run without any input if none is specified, defaulting to an IBM with a population size of 100 for one generation. All other options, including non-standard functions, are left to the user.

Working this through, I’m slightly apprehensive about the motivation for including gmse_apply function. Once you strip the mechanistic approach from the resource and observation models, all you really have are two values: (1) the population abundance or density and (2) the estimate of the population size or density. Once you include the manage_target into gmse_apply (necessary, I believe), then the genetic algorithm is really just a fancy way of getting the difference between the population estimate and the target size, and then setting a number of culling actions acceptable for users. Users then cull as much as possible because they’re assumed to want to use the resource as much as possible. Of course, we can consider other parameters that affect user actions (e.g., maximum catch, effort), but if we’re interested in learning about how these concepts affect harvesting in theory, then they can and should probably be studied using a simpler model. The real point of the genetic algorithm is that it allows for complex, multi-variable goal driven behaviour, as might occur given indirect effects (e.g., organisms on crop yield) or multiple options (e.g., culling versus scaring or growing) and spatial complexities. There seems little to be gained by calling the genetic algorithm to tell users to cull as much as possible, which can be done with a (very) simple function.

Update: 28 AUG 2017

I have finally fixed the annoyance in the shiny app of GMSE that caused the bottom of the browser to black, hence making it difficult to set parameter values in some tabs.

Additionally, by hovering over the different options in the application, the software user can now see a brief description of what each option does in the simulation.

Update: 23 AUG 2017

I am experimenting with ways of demonstrating the conflict between what a manager incentivises, and what the users actually do, in GMSE. Below are some plots that show this for a few sample simulations. The five panels in each plot correspond to the five possible actions where policy is set. Policy set by the manager is shown with the black solid line, with the thin coloured lines reflecting individual user effort expended into each action.

The right axis is fairly easy to interpret – it’s just the percentage of the user’s total budget devoted to a particular action (note, this is not necessarily the number of actions a user performs because different actions can cost different amounts – hence the term ‘effort’).

The left axis is a bit trickier – it’s how permissive of an action the manager is in practice. High values correspond to an action being highly permitted by the manager (i.e., the manager invests no effort in making these actions costly), whereas low values correspond to an action being less permitted (i.e., the manger invests highly in making these actions costly for users).

The end result is that the lines indicating manger permissiveness are typically correlated with user effort towards any particular action. In the first example below, this is true for scaring and culling (as the manager becomes more permissive of these actions, users tend to take advantage and spend more effort doing them). Note that users do not feed because they have nothing to gain by feeding the resources, even though the manager is usually permits feeding (around generation 75, the population started going way over the manager’s target).