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:

Update: 03 AUG 2021

A new version of GMSE 0.7.0.0 has been released on CRAN now. The changelog provides a summary of the changes and bug fixes. One new feature was written by Adrian Bach, and now GMSE includes the ability to use trajectory-based prediction to determine whether or not intervention on the part of the manager. The new feature will be used in a manuscript in preparation (note that on the GMSE website, I have now included a list of publications that use GMSE). This includes three new arguments:

• traj_pred is an argument that takes a TRUE or FALSE value that determines whether or not the trajectory of a population should be used when a manager is deciding to intervene. This argument applies only when action_thres > 0 (i.e., when there is a range around the target population density within which a manager will not intervene by updating policy). When traj_pred = TRUE (default is FALSE), the manager will use a simple linear extrapolation from the density of the population in the previous time step and the current population density to predict what the population density will be in the next time step. If the predict density is outside of the acceptable range specified by action_thres, then the manager will intervene.

• bgt_bonus_rest is an argument that takes a TRUE or FALSE value and determines whether (TRUE) or not (FALSE) the manager’s budget bonus for inaction will be reset to zero following a time step of policy update. If FALSE, then budget is only reset when costs of culling has been decreased in the previous time step. We will probably want to generalise this to the cost of all actions later, which will involve doing a bit of work to this block of code. Else we might consider renaming the argument to indicating culling somehow?

• mem_prv_observ This parameter is an argument that takes a TRUE or FALSE value. When TRUE, it allows the manager to remember the population size from the previous time step. It should be used with the traj_pred option.

Regarding the last bullet, I’ve made an update to GitHub v0.7.0.1, which fixes a minor bug that I just noticed this morning that unfortunately made its way into the CRAN version. I don’t think that it will actually affect any simulations though, so I don’t see a need to immediately update on CRAN. See Issue #76 for more details.

There are new options introduced by Jeroen for incrementing user and manager budgets using yield from the landscape.

• usr_yld_budget takes a numeric value that acts as a multiplier on the user’s total yield. The value of the argument is multiplied by the user’s total yield to get the user’s budget increment. If, for example, we wanted to simulate a situation in which users’ budgets were entirely based on their yield, we could set user_budget = 0 and user_yld_budget = 1. Note that regardless of yield, user budget can never be less than 1 or greater than 100000.

• man_yld_budget is the equivalent to usr_yld_budget, but for the manager. Manager’s yield is accumulated on public land, so if public_land = 0, then there will be no increment to the manager’s budget.

Improvements to the gmse_gui have also been made by Jeroen, which can be accessed online here. Note that I have not yet included the trajectory strategy options that Adrian developed into the GUI.

After some conversation, I decided to change the movement operations when resources feed more than once (times_feeding > 1) or require a threshold amount of food for survival (consume_surv > 0) or reproduction (consume_repr > 0). Previously, resources would feed, then move times_feeding times. After this moving and feeding, resources would move again at the end of the time step. Somewhat confusingly and arbitrarily, resources would therefore get two moves in a row under these conditions before the timestep ends. In GMSE v0.7, resources now only move once after feeding.

The new GMSE version also fixes some minor bugs, most notably getting rid of the rare cases in which resources were able to feed once on the landscape after dying. There are also new error catches for gmse and gmse_apply functions to make sure that arguments are structured correctly (e.g., quickly flags when a list or vector is used instead of a scalar value).

Update: 21 MAY 2020

Several new features have been added, which will also make their way into the next version of GMSE v0.6.0.x. Two new arguments usr_yld_budget and man_yld_budget are now available in gmse and gmse_apply. These arguments allow an increment of user and manager budget from yield obtained from the landscape. Total yield from one time step of for a user is multiplied by usr_yld_budget, and this product is then added to whatever the user’s baseline budget was. Similarly, the mean total yield of all users is multiplied by man_yld_budget, and this product is then added to whatever the manager’s baseline budget was. There was a bit of an issue when integrating the code with the action threshold budget bonus, so the solution was to add two new columns to the agent array. The budget for any agent is then just the sum of the basline and each of these columns. Hence future development affecting the budget bonus and yield to budget effects can be more easily done independently.

A new vignette has also be written by Adrian, which will be on the updated GMSE website soon.

Update: 18 MAY 2020

A bug has been found and fixed, which was introduced into gmse_apply through the changes to the perceptions of actions for users and managers in commit 86c298eddd1c778defaac04d30ab342eef58fb50. The problem was that the internal elements of the paras vector were changed to reflect columns where perceptions of different actions were to be located in gmse. But in gmse_apply, the same change was not made, so that all of these elements were interpreted as low integers; it’s remarkable that a segfault did not occur in any of the tests. Jeroen pointed out a major difference between the dynamics in runs of gmse and gmse_apply given identical parameter values, which led the investigation into what was causing the problem. Finally, we found the aforementioned problem, which was resolved in commit 7a97eb08171147f688036851be5e3783ab631495. A subsequent fix was made to the testthat function paras vectors where needed. The dynamics in gmse and gmse_apply are once again the same. On Friday, we have tentative plans for a crash party with ConFooBio members to download the newest version of GMSE on the rev branch and try to crash it.

Update: 17 MAY 2020

Over the last two days, I have made some changes to the code and done some additional testing. Changes are summarised below.

Fixed memory leak

Since I had been doing some substantial development in C, I thought it would be wise to run valgrind to check for memory leaks. I found one, due to my forgetting to free some memory in the new landscape.c file. The memory is now correctly freed in the file, and running gmse and gmse_apply in valgrind produces no memory leaks.

==33279==
==33279== HEAP SUMMARY:
==33279==     in use at exit: 67,085,047 bytes in 12,107 blocks
==33279==   total heap usage: 896,509 allocs, 884,402 frees, 557,718,069 bytes allocated
==33279==
==33279== LEAK SUMMARY:
==33279==    definitely lost: 0 bytes in 0 blocks
==33279==    indirectly lost: 0 bytes in 0 blocks
==33279==      possibly lost: 0 bytes in 0 blocks
==33279==    still reachable: 67,085,047 bytes in 12,107 blocks
==33279==                       of which reachable via heuristic:
==33279==                         newarray           : 4,264 bytes in 1 blocks
==33279==         suppressed: 0 bytes in 0 blocks
==33279== Reachable blocks (those to which a pointer was found) are not shown.
==33279== To see them, rerun with: --leak-check=full --show-leak-kinds=all
==33279==
==33279== For counts of detected and suppressed errors, rerun with: -v
==33279== ERROR SUMMARY: 2736 errors from 4 contexts (suppressed: 0 from 0)

I have also initialised a new file in the notebook for functions that we will want to periodically run through valgrind.

Perceptions of actions for users and managers

Following a conversation with Jeroen and Nils, we decided that another change would be beneficial for GMSE v0.6.0.x. In all previous versions of GMSE, how managers and users perceive the consequences of scaring actions is hidden. The consequences were calculated in a reasonable way by GMSE and used within the fitness functions of the genetic algorithm to determine which strategies of setting policy (for manager) or acting (for users) was best. The relevant values for the perceived effects of scaring, culling, castration, feeding, helping offspring, tending crops, and killing crops were set partly in the paras vector, and partly calculated internally in C. Now, instead, they are calculated for each agent (manager and users) and placed directly in the AGENTS array.

The AGENTS array has therefore increased by seven columns. Rows in each of these columns hold the perceived effects of each action for each agent. I’m not sure why I did not think of this before, but there are several benefits to doing it this way. First, anyone running gmse or gmse_apply can now directly alter how users perceive the effects of their own actions using the perceive_scare, perceive_cull, perceive_cast, perceive_feed, perceive_help, perceive_tend, perceive_kill arguments, which set how perceived utility will change given a single instance of scaring, culling, castration, feeding, helping offspring, tending crops, and killing crops, respectively. By default, all of these arguments will be calculated in a realistic way, as has been described already for the developing version, so the arguments can be ignored and GMSE will run as normal. But if, e.g., there is reason to model the users’ perceived efficacy of culling as being much higher than that of scaring (e.g., based on parameterisation from empirical social data), then more appropriate values for a specific system can be used. Additionally, with gmse_apply, the agents array can be changed directly, meaning that manager perception of action effects can be changed, and that perceptions among users can potentially be varied (e.g., if some users are more convinced that scaring will work, while others are more sceptical). I also like that the perceived effects are being returned in the agents array, so the whole parameterisation of the fitness function is much more transparent now.

The code now runs with no memory leaks, and GMSE builds on the rev branch with no notes, warnings, or errors.

Update: 15 MAY 2020

Minor bug fix

Jeroen and I have fixed a minor bug in the GMSE code that was causing a crash given the very restricted situation in which the package is built using the Rstudio drop-down menu. In such builds (but not, curiously, when building from the console), a fatal crash can be caused when looping the gmse_apply function with variable names that are identical to argument names (e.g., stakeholders = 8, followed by gmse_apply(stakeholders = stakeholders)). Just in case, to avoid this issue entirely, I have inserted a short function into gmse_apply that produces an error when attempting to set variable names equal to argument names. This can be completely circumvented in gmse_apply by setting the hidden argument my_way_or_the_highway equal to anything (i.e., not NULL). I have also included a garbage collect at the end of gmse and gmse_apply functions.

Consistency between gmse and gmse_apply

In previous versions of GMSE, the functions gmse and gmse_apply had a few slightly different defaults. This sometimes led to slightly different dynamics when running simulations and comparing between gmse and looped gmse_apply under default arguments. To avoid any cofusion, Jeroen and I have decided to make gmse and gmse_apply default arguments completely identical. The equivalence of the two functions given the same initial conditions has been confirmed through simulations.

Update: 12 MAY 2020

Updated vignettes

I have now checked vignettes 5-7 to see if updates were needed. Only an update to vignette 7 was necessary. Note that even where the code for vignettes remains unchanged, it would still be a good idea to re-knit everything for the website as a test to make sure that everything runs as expected with new defaults for v0.6.0.x.

Uneven land allocation

Jeroen has incorporated some code and a new ownership_var argument to adjust the variation in land ownership among users. I’ve re-worked the way that public land is handled in the new landscape owner building function to account for this, and everything appears to work as intended.

Update: 11 MAY 2020

Bug fix to public land allocation

A small bug was noticed by Jeroen this morning relating to public land allocation. Essentially, when the proportion of public land amounted to fewer cells on the landscape than allocated to a single user, public land would default to zero cells. We needed some way for there to be nonzero public land that was also quite small. Given such situations, a new function small_public_land is called. Public land is placed in the centre of the landscape, thereby cutting into at least two users’ owned cells equally. First a rectangle proportional to the landscape is produced, followed by more public land cells attached to each side of this rectangle as necessary to achieve approximately the correct number of cells. It is important to recognise that simulation parameter values need to be chosen with care. If the total number of landscape cells cannot be evenly distributed among users and public land, then of course deviations from gmse and gmse_apply arguments will occur. The updates cause small deviations to occur for larger sections of landscape when possible. For example, the below produces a landscape with four equal quadrants, three for users and one for public land.

sim <- gmse(time_max = 4, stakeholders = 3, public_land = 0.25,
land_dim_1 = 100, land_dim_2 = 100, land_ownership = TRUE);

The below is a situation in which the amount of public land is not small. A total of 2000 public land cells are called for, with the remaining 8000 to be divided evenly amongst nine users. This is impossible, and the shortest-splitine algorithm cannot give that level of precision and make the landscapes clean.

sim <- gmse(time_max = 4, stakeholders = 9, public_land = 0.20,
land_dim_1 = 100, land_dim_2 = 100, land_ownership = TRUE);

Instead, 1840 cells are allocated to public land, with each user being allocated rectangular blocks of 900 or 920 landscape cells. But when public land to be allocated is very low, fewer cells than given to a single user, deviations are then placed on users to make the public land proportion as accurate as possible in the centre of the landscape.

sim <- gmse(time_max = 4, stakeholders = 7, public_land = 0.03,
land_dim_1 = 100, land_dim_2 = 100, land_ownership = TRUE);

In the case of the above, e.g., exactly 300 public land cells are put down in the centre of the landscape, while users are given somewhere between 1303 and 1462 cells.

Update gmse GUI

The gmse_gui function has now been updated so that all options are available (and default landscape allocation is enacted) for running simulations from within the shiny GUI. I have not put this up on the shiny server yet, but on the rev branch, the new GUI can be used within R. This now concludes all of the planned updates for code changes to v0.6.0.x, and now the last two things to do will be to add vignettes and update the website. Please contact me to request any additional changes.

Update: 10 MAY 2020

Yesterday and today I implemented the fourth bullet point listed on 6 MAY 2020, which changes the default pattern of land ownership in GMSE. In over a decade of coding individual-based models, I think that this might be the first time that I have needed to apply recursion to a problem. I actually only just learned that recusion is possible in R, but decided to write the function in C just to get a bit more speed.

The goal of the new algorithm is to divide a rectangular landscape of length land_dim_1 and width land_dim_2 into stakeholder smaller rectangles that are as equal as possible in area (i.e., cell number). To do this, I used a shortest-splitline algorithm in a new file landscape.c. The relevant function is the break_land function. For a given rectangular landscape, the break_land function draws a line splitting the landscape along its longest dimension (i.e., if the rectangle has a length of 20 and a width of 10, it is split lengthwise). If the number of owners to be placed within the landscape is even, then the line splits exactly halfway down the middle of the landscaope. If the number of owners to be placed within the landscape is odd, then the line is placed in proportion to the ratio of ownership – in other words, if 3 owners need to be placed on the landscape, then the split line would be placed 2/3 of the way along the longest dimension (3/5 of the way for 5 owners, 4/7 of the way for 7 owners, and so forth). The rectangular landscape is then split into two smaller sub-rectangles. The same splitting procedure is then run on these sub-rectangles with the new number of owners (i.e., half of the original in the case of an even split, and appropriate portions such as 1 and 2, or 2 and 3, for the odd number of original owners, with the biggest rectangle getting the higher number of owners). The splitting of the landscape into smaller and smaller sub-rectangles continues until the number of owners that need to be placed on a rectangle equals one, in which case an integer representing an agent’s ID is placed on all landscape cells in the rectangle.

In the case of public land, the correct proportion public_land is put on the landscape by running the above algorithm in the break_land function with an appropriately added number of owners. If, for example, half of the land is meant to be public, then the total number of unique sub-rectangles to be made is adjusted to users / (1 - proportion_public_land). Where if proportion_public_land = 0.5, e.g., the total number of sub-rectangles to be made will be twice the number of users. Once all of the cells are assigned an integer value by the algorithm, all cells with a number higher than any user ID are then re-assigned as public land.

These changes have now been pushed to the rev branch, so as of writing this, they are available to use. Relevant testthat functions have been adjusted where necessary to account for the code changes, and the new code passes all checks and all CRAN tests with no notes, warnings, or errors.

I like the above algorithm so much that I do not see any real need to even make the old straight-line simple sort approach an option in gmse. The old approach was a simplification that always seemed a bit unrealistic, and I always wanted to come up with a better way of allocating owned landscape cells. Please contact me or raise an issue on GitHub if you really want the old way back in gmse (note that gmse_apply of course always allows custom landscape definition).

With this new feature, I have now tackled bullet points 1, 3, and 4 from 6 MAY 2020. The next thing that I will do will be to update the GUI, which I do not expect to take long. Once this is finished, I will begin writing my vignette on resource density, and updating the website accordingly. I think I might also update the vignette on default GMSE structures to briefly discuss the above algorithm and how landscape ownership is assigned. I will put out another announcement before sending to CRAN, but please let me know if there are changes that you would like to see in GMSE v0.6.0.x.

Update: 8 MAY 2020

I have made some progress on the plans from 6 MAY 2020. None of my updates are incorporated into the main R branch, so they are not permanent and I am still happy to take suggestions for improvements or other ideas in addition or in contrast to what folows.

I have completed the first bullet point of the plans from 6 MAY 2020. The order of operations is therefore changed in the resource model on the work-in-progress rev branch. Hence, calls to the default resource model proceed by first having resources consume on the landscape, then birth, death, and movement. Within resource consumption rules, feeding happens first, followed by movement for up to times_feed iterations. The landscape consumption also now calls either one function or another: one function where consume_surv > 0 or consume_repr > 0 in which consumption happens in a random order among resources for a variabl enumber of times_feed, and a second (faster) function where consumption happens once for all resources (in order of the array, as the order of consumption does not affect the resources themselves).

Despite the information in the vignette on GMSE structures, column 13 of the agents array was actually being used within the more or less unused anectdotal function, in which each agent ‘looks around’ from the cell on which they are standing to see how many resources are in the general area. This function has no purpose at the moment (I was thinking about using it to model disagreement with the manager’s estimate and the users’ inferences about population size), but writes to the column anywa in gmse (but not gmse_apply). I’m going to refer to it as written to, but unused in the vignette, and use column 14 of the agent array to record the number of cells owned.

Writing to the number of cells owned is done in a new count_agent_cells in R or the count_owned_cells function in C. I wrote the R function first, but then had second thoughts about relying on it in gmse and particularly gmse_apply. What if we later want land ownership to change over the course of a simulation? I’m not sure in this case, actually, how much extra speed we get from C (it relies on two nested for loops across the landscape, rather than a which or == in R, which I suspect is actually just as fast since the underlying C is probably better than mine). In any case, it seems to make sense to have the owned cells calculated at the beginning of the user function, so this is what I ultimately decided on. But the R function can also be used.

I have changed the fitness function for users as such:

• Scaring is now assumed to decrease the number of resources affecting the user by one minus the proportion of landscape cells owned. Hence, when less land is owned, scaring is assumed to have a bigger effect because resources will more likely move to someone else’s land. Also, if no land is owned, scaring is assumed to have no effect, for obvious reasons.
• Culling is now assumed to decrease the number of resources by one plus the expected number of offspring the resource is assumed to produce in the time step. This expected offspring number E_off is the lambda parameter, plus the amount a resource consumes on a landscape cell divided by the amount needed to be consumed for producing one offspring.
• Castration is assumed to decrease the number of resources by E_off.
• Feeding is assumed to increase the number of resources by E_off.
• Helping offspring is assumed to increase the number of resources by one.

These changes are now pushed to the rev branch, which is still under ongoing development. This effectively takes care of the third bullet point from 6 MAY 2020. Hence, bullets one and three are now done, so here is what remains.

• Update the GUI
• Add new landscape patterns of land ownership options
• Two new vignettes to add to the package and website
• Other website updates that need to change (I will actually need to re-read everything so that any old information that is now wrong is replaced from existing vignettes and documentation.

I will try to make the above changes as soon as possible, but the changes that will have the biggest effect on simulations are now finished, so people can play around with the new options. The update in progress builds with no warnings, errors, or notes, and passes all tests in testthat.

Update: 6 MAY 2020

New version of GMSE

A new version of GMSE v0.5.0.0 is now available on GitHub. Here are the major updates:

• As mentioned on 23 APR 2020, new arguments for implementing manager action thresholds and budget bonuses (for not acting in a time step) are available using the arguments action_thres and budget_bonus in gmse and gmse_apply.
• Also as mentioned on 23 APR 2020, a minimum age of resource reproduction can now be set using the age_repr arguments, and random sampling from a range of user budgets around user_budget can be set to be plus or minus usr_budget_rng in gmse and gmse_apply.
• More complex interactions between resources and landscape can be implemented, such that resource birth and death depend on resource consumption. These values are set using the new arguments consume_surv, consume_repr, and times_feeding. See details below.

The new version can be downloaded from GitHub, but will not be submitted to CRAN. It is likely to be the last version of GMSE that is completely backwards compatible with previous versions. In other words, if you set parameter values for previous versions of GMSE as arguments in gmse or gmse_apply and ran simulations, then you are effectively running the same code (and should get the same results) in v0.5.0.0. In the next version update, this will not necessarily be true because minor changes to the default order of operations are planned (see below). In practice, this is highly unlikely to change any results from simulations in previous versions, but it is good to be aware of what has changed if you are expecting old code to run identically on new versions of GMSE.

Plans for GMSE v0.6.0.x: what will change

While backwards compatibility has been a priority up until now with GMSE development, the most recent developments in GMSE v0.5.0.0 make it clear that the benefits of make some small changes to the default simulations outweigh the costs. This mostly has to do with the order of operations in the default resource submodel (currently: movement, feeding, birth, death). Given the new option of resource consumption affecting birth and death rates, and the consequences of farmer actions (i.e., when land_ownership == TRUE), a change in the order makes sense. When feeding affects birth and death in v0.5.0.0, resources first move once, then move again, then eat, move, eat, move, etc. for as many times_feeding as defined. Hence, resources are moving twice at the start of a time step, and potentially landing on good landscape cells, but then not eating at all. With respect to farmer actions, if resources move far in a single time step, then their efforts might have very little actual effect on their individual yield. This is because, since movement happens first, resources might be highly likely to move to another farmer’s land anyway, making the actions futile. This might not be a problem for modelling some systems (where such actions are in fact futile), but the fitness function of the genetic algorithm is designed in such a way as to have farmers assume that one scare or cull will decrease the effects of a resource accordingly. In practice, this doesn’t matter too much because the absolute values in the fitness functions are arbitrary for now (i.e., we could multiply them by the probability that a resource did not leave the farmer’s landscape and expect the same farmer actions – there’s nothing else they can do). But it could become an issue later, and we might at some point be really interested in looking at patterns among different users in a more interesting way. Therefore, in GMSE v0.6.0.x, I plan to change the order of operations in the default resource submodel (some more of my reasoning is discussed in Issue 56). Overall, I plan for the following changes:

• Change the order of operations in the resource submodel to resource feeding, birth, death, movement. In feeding subfunction, the order will also be changed to first feeding, then movement. Hence, the first thing that a resource will do at the start of the time step is feed, then move, then feed, etc., until they have exhausted their times_feeding. If the resource is not scared or culled at the end of the previous time step (user submodel), then it will feed on the landscape cell on which it last resided; scaring and culling will therefore directly affect a farmer’s landscape. At the end of the resource submodel (after all of the feeding, birth, and death have happened), resources will once again move. This is because the feeding process is not forced in GMSE, so if default consume_surv = 0 and consume_repr = 0 were used, movement would not happen at all in the resource submodel. These changes will be small in the code, but will change the default behaviour of gmse() and gmse_apply().
• Update the options for the GUI when running gmse_gui(). This is just an oversight from the last version, as I should have done this for the new arguments available in v0.5.0.0. The new GUI will look the same as the old (I’m not sure if anyone is actually using it, but it’s worth updating just in case), and I will send it to the same place on shiny so that it can also be run without using R.
• Change the fitness function for users here in the genetic algorithm. Users assume that scaring or culling has an effect of removing one resource, that castrating has an effect of decrementing by lambda, feeding increments by lambda, and helping offspring increments by one. This has worked for a while, but I think some modifications will make it more realistic.
• Scaring should have a predicted effect of decreasing a resources by 1 - land_owned / total_land. If one farmer owns all of the land, then scaring doesn’t make any sense, and this will be reflected in the user fitness function. Similarly, when land_ownership = FALSE, scaring also won’t make any sense (previously, people using GMSE just had to recognise that the combination land_ownership = FALSE and scaring = TRUE was nonsensical, now the users will recognise it themselves).
• Culling should really have an effect of 1 + E[offspring]. That is, 1 + lambda where consume_repr = 0, and something else where consume_repr > 0. In the absence of any better ideas, I propose that the user assumes that a resource will eat all of the contents of a single landscape cell and produce the consequent number of offspring. This will be false most of the time, but it seems reasonable from the perspective of a landowner, who would have no way of knowing how much more a single resource would be expected to eat and reproduce as a consequence of feeding. Hence, from the user’s perspective, the effect of culling will be a decrease in resources of 1 + lambda + (res_consume / consume_repr).
• Castration, along the same lines of culling above, should decrease resources by a magnitude of lambda + (res_consume / consume_repr).
• Feeding should, similarly, increase resources by a magnitude of lambda + (res_consume / consume_repr).
• Helping offspring can be kept as incrementing resources by one.
• Default patterns of land ownership on the landscape can be changed (as long as we’re making changes that affect backwards compatibility). The long thin bands of land have worked up until now, but, especially when the number of land owners is medium sized relative to the number of cells (e.g., 100 users on a 100 by 100 landscape, giving each user a 100 by 1 strip of land), this doesn’t look quite right. Instead, a new R function can make the default option a polygon design using a shortest-splitline algorithm. This will look a bit nicer and more realistic.
• Two new vignettes, the first written by Adrian explaining the action threshold and budget bonus options. The second written by me explaining resource density regulation (including the new consume_surv and consume_repr arguments).
• Any other updates to the website that will inevitably need to be made.

The end result of all of this will be, I think, an improved GMSE, but also one that breaks from previous versions in the default options and potentially stakeholder behaviour. For example, land owning users will inevitably see culling as a slightly better option for increasing yield, all else being equal (i.e., given identical costs, as set by the manager). I’m keen to start this new version as soon as possible. The new version will be pushed to CRAN after these changes have been implemented. I will attempt this by the end of the month, or early June at the latest. Please contact me either via email or GitHub if you have any thoughts or concerns about the next version.

Update: 23 APR 2020

With commit 6d1e5e6ca41bc6bb45c62cb174f34480c3c5add1, I believe that I now have resolved Issue #56 and created a version of GMSE that now allows resource death and reproduction to be tied to resource consumption on the landscape. Initial testing shows that this works in both gmse and gmse_apply. I have not had a chance to code in the suggestion from Jeroen to have it be the probability of survival or reproduction that changes with increasing consumption. Rather, given that a resource has consumed a sufficient quantity of resources consume_surv, the resource will survive. A separate new gmse argument consume_repr ties reproduction to landscape consumption by having resources produce floor(consumed / consume_repr) offspring given some total consumption consumed.

There is still considerable stochasticity introduced in terms of whether or not a resource will survive or reproduce. First, each resource can potentially feed more than once by using the new times_feeding argument. When times_feeding > 1 (default is 1), resources will consume the landscape, then move, then consume, and so forth until all of their times_feeding have run out. In practice, all resources get the same times_feeding value, but the code for executing it depends on a new column in the resources array. This means that there can be variation in the number of times a resource feeds if using gmse_apply (or gmse, if we later decide to sample feeding rates from a distribution). Feeding events also occur in a random order, which avoids the first resource on the list filling themselves up before moving on to the next resource (as would be a natural way to program it with a for loop). I did it this (slightly less efficient) way to more realistically model how resources would really move around simultaneously on a landscape eating as they move from place to place.

Hence, in the end, there is a lot of uncertainty surrounding whether an individual resource will have enough to survive and reproduce in a time step. But the probabilistic idea of Jeroen could also be imposed with a function of some sort that would link consumption to a probability of death (will need to think of the best way to actually do this in the code though). Testing shows that resource increases to a carrying capacity imposed by the constraints of the landscape, with a lot more stochasticity than occurs just by settting a value for global res_death_K. To get started, for example, the below demonstrates the new features.

sim <- gmse(consume_surv = 2, consume_repr = 3, times_feeding = 6,
res_birth_type = 0, res_death_type = 0, land_ownership = TRUE,
stakeholders = 12, scaring = TRUE, time_max = 40, user_budget = 1);

I set the user_budget to unity above just to show how the population rises to a natural carrying capacity when the manager and users are unable to cause any population change. We can recover a (usually) well-regulate population by settting the user_budget back to the default.

sim <- gmse(consume_surv = 2, consume_repr = 3, times_feeding = 6,
res_birth_type = 0, res_death_type = 0, land_ownership = TRUE,
stakeholders = 12, scaring = TRUE, time_max = 40,
user_budget = 1000);

Just to clarify, each landscape cell still puts out a value of 1. The argument consume_surv = 2 means that resources need to consume two cells worth to survive a time step (by default, resources consume half of what’s on a landscape cell with each visit). The argument consume_repr = 3 means that resources need to consume three cells worth to produce one offspring. The argument times_feeding = 6 means that resources get to feed on six cells in each time step.

Setting res_birth_type = 0 and res_death_type = 0 effectively removes all other factors in GMSE that affect resource birth and death. Hence, what is observed from running the above code is purely a consequence of the new rules limiting survival and reproduction due to landscape consumption. Note that these values need to be specified, else GMSE will assume that other defaults apply in addition to the effects of landscape consumption. In other words, the global external carrying capacity from res_death_K will still limit population growth because by default res_death_type = 2, and the mean birth number lambda will still apply because res_birth_type = 2. I don’t want to override this just in case anyone wants to include multiple types of processes affecting birth and death (I also want the new version to be backwards compatible with existing code), but I think it will be important for me to write some documentation that makes using these new features really clear. I’ll start a new vignette.

The new version is GMSE v0.5.0.0 with Adrian now included in the list of contributors. The action_thres and budget_bonus can now be run in gmse and gmse_apply, and usr_budget_rng and age_repr are also available options. The new version passes all CRAN checks with no notes, warnings, or errors. The new version should be completely backwards compatible. If you have existing code to run, it should do the same thing in the new version that it did in previous versions. Everything is now just a new feature that can be used, but is not by default.

The last thing that I think I want to do before merging with master, updating the website, and sending to CRAN, is to make a new default landscape. It would be good to make each landscape into more farm-like blocks rather than the narrow stripes that exist now. I don’t have a good algorithm for doing this for an arbitrary number of stakeholders (i.e., for N stakeholders, separate an X by Y landscape into roughly equal blocks of sizes ca (X*Y / N)). If anyone wants to have a go at this and create a function that does it, please let me know! Else I’ll try to get around to it before early next month where the plan is to send to CRAN. If anyone has time (Nils Jeroen Adrian), feel free to download from the dev branch and go.

install.packages("devtools");
library(devtools);
install_github("ConFooBio/GMSE", ref = "dev")

It should be up now.

Update: 13 APR 2020

I have now completed two major tasks in a GMSE update, first mentioned earlier in the month. The new code is up on the rev branch.

I have resolved Issue #53, so partial matching is not allowed in gmse_apply (or in the tests). I have also resolved Issue #59, meaning that it is now possible to use different action thresholds and budget bonuses in both gmse and gmse_apply (arguments action_thres and budget_bonus), and variable user budgets can now be specified in gmse, and in gmse_apply (argument usr_budget_rng – of course, manual edits to user budgets were always allowed in gmse_apply by directly editing the agents array).

In making these changes, I have also allowed gmse_apply to keep track of the current time step when old calls to the function are passed back through old_list (i.e., when looping old list as in here, the time step will also be tracked automatically through the loops, not just outside of gmse_apply through the loop index – this makes the action threshold possible, as it only applies after the first time step).

I now need to move on to Issue #56, which will more or less complete the new version, though I also might update the default landscape ownership placement.

New tests have been written in testthat, and the package builds with now errors, warnings, or notes.

Update: 1 APR 2020

I am now beginning a moderate update to GMSE, which I hope to have complete and up on CRAN by the end of the month. There are three specific goals of this update:

1. Resolve Issue #53: Just to head off any potential issues in gmse_apply, several places in the code use the format sim_list$element_name rather than sim_list[[“element_name”]]. The latter is more robust because it prohibits partial matching. Hence in gmse_apply.R, this should be changed. 2. Resolve Issue #59: The code from AdrianBach needs to be incorporated into the master branch of GMSE. Doing this will mainly consist of resolving a merge conflict and a bit of refactoring, and I plan to do this early on in developing what I hope to be a new version of GMSE by the end of APR 2020. 3. Resolve Issue #56: This is the major one. I want to add options in GMSE to allow resource dynamics to be dependent up on landscape properties. Specifically, for resources to consume the landscape cell contents, and for options to exist that allow resource birth, death, and (potentially) movement to be dependent on consumption. This will mainly involve updates to the resource model, with updates to the gmse and gmse_apply function as needed. 4. In addressing 1-3, I will need to check to make sure that gmse_apply is handling all of the new options correctly, and new tests will need to be written as necessary. I will continue to update as need be and make sure that any updates from other team members (e.g., Issue #58) are merged appropriately before the new GMSE version. Update: 4 FEB 2020 Yesterday, I received an email noting a problem with the GMSE documentation. Specifically, a warning arose because ... was in the documentation but is no longer allowed as an argument in the default R functions. Notice the issue on this line of the manager.R file, which is contradicted later by the function itself. This has now been resolved. Update: 15 MAR 2019 GMSE v0.4.0.11 is now available on GitHub. Some minor changes and bug fixes are included in the update version. • Issue #46 appears to be resolved. This was caused by some memory management problems, making it impossible to increase the paras vector in development. It also made repeated calls to gmse_apply crash when the exact same arguments were specified but some values (e.g., resource abundance) changed. • A new section on the website and in the corresponding vignette explains how to loop custom sub-models in gmse_apply, noting some special considerations that could potentially cause confusion. • As a consequence of the change in the way that base R generates random numbers, tests using the testthat needed to be re-worked. This has been done, and GMSE v0.4.0.11 should be on its way to CRAN as of next week. I will update once more when GMSE v0.4.0.11 is successfully submitted to CRAN. Update: 21 AUG 2018 GMSE v0.4.0.7 is now available on GitHub, and the official GMSE repository is now transferred to the ConFooBio organisation. I am currently in the process of fixing some website issues, but most of the important stuff has transferred to the new website location. 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, if 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 too 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 by 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 only 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]; adjusr <- manager_acts[estim_row, 5]; observations[i,j+1] <- target - adjusr; #---- 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).