confoobio / gmse Goto Github PK

Just allow an argument in gmse_apply that allows the full results of a previous gmes_apply to be passed back into it. Then force an error when someone tries to pass a PARAS.

gmse() crashes R when either the manager and user budget is zero. Also, cannot have cost values of zero (fixes in initliase).

For some reason, having more stakeholders appears to lead to less culling of resources even when all of them are attempting to do it. If there are more stakeholders to act, then actions should happen more often because each has the same budget.

Check to ensure that the ownership of cells on the landscape is being divided accurately -- at the moment, ownership is being assigned probabilistically, not exactly (e.g., expected that two owners each get half of the land, but some variation around this expectation). It might be a useful option to allow exact specification of land ownership cells.

For some reason, the user.c call appears to crash when there is more than one resource. I only noticed this after the overhaul of the observation model, but I doubt they are related. More likely, one of the many arrays with dimensions that depend on resource number is being built or called improperly, leading to a segmentation fault. This should be fixed, of course, but I suspect the problem is not too far buried.

Currently, only the history of interactions from the previous time step directly affects the genetic algorithm for stakeholders and managers. For managers especially, this could be made a bit more nuanced. The entire history of total actions and resource dynamics is recorded, and this could easily be made available (e.g., in PARAS_REC) for managers to make decisions. Incorporating these data into the genetic algorithm, and therefore into agent decision making, could be tricky, but one simple example of this could be having managers use the per-time step mean number of stakeholder actions in the last 2-3 time steps to predict future user actions with a bit more inertia. Managers could also use stakeholder action history from earlier time steps, but weighting each by how long ago they occurred.

Currently, to simulate mark-recapture observation and analysis, values for fixed_mark and fixed_recapt need to be specified in GMSE, and the manager would have exactly these numbers of marks and recaptures in each generation, respectively. It would also be useful to, instead of specifying exact numbers, to have the manager search a general area, then mark all resources in that area. Next, the manager could search again and recapture, so the exact number is not always set and the observation process probably mimics more closely what happens in the field. This type of sampling is actually already available (observe_type = 0), so I would just need to add some code to have managers interpret some observations as marks and others as recaptures.

In case two of the observational model, the length and width of a sampling block both increase linearly with the number of agents doing the sampling; hence, sampling area increases exponentially with the number of observers, which is probably unrealistic. There are two ways to potentially address this:

1. Simply add another square for each agent observing, in the next place it would otherwise go.
2. Probably more easy, have a countdown defined by += (int) agent_array[agent][8]. Only allow resources to move when this countdown hits zero, and reset it it thereafter. Hence, observers will observe more n more blocks if there are n more observers.
3. Could adjust the squared dimensions appropriately to retain a square block, but this seems inefficient, imprecise (due to rounding), and unnecessary.

Related to discussion with Jeremy regarding the Islay geese, need to have a kind of observational model in which agents move to take measurements, but resources move along roughly the same time scale. This can of course be accomplished one way if we:

• add an option for resource movement during the course of observation: Do this in the main observation function as an if(resource_movement == 1) type criteria at the tail end before the break (to avoid unnecessary movement). This will require also including the resource movement function (currently in resource.c) in the observaction.c file. May as well just dump the whole thing in in the interest of modularity, though if it stays the same, it will be tempting to create a utils.c file of some sort. This resource movement option can be applied to the existing method case 0, as appears in the switch function of the observation function.

To do a sweep of the landscape while allowing resources to move, I think we'll want a completely different method of population size estimation. What this method will do is:

1. Start an observer(s) in a fixed x location of x = 0 on the landscape
2. Have the observer(s) census all individuals on the x locations x to x+view (i.e., observe view rows)
3. Set a new x = x+view
4. Move resources
5. Iterate 2-4 until x+view is greater than the y dimension land_y
6. If there are any ys left, then iterate the last x to land_y.

The procedure above will simulate observations over a time that is proportional to their view (and thus ability to census) -- the more time it takes, the more the resources can move and potentially lead to measurement error. The observational array returned will still be output in the same way -- resources will be marked as with the case 0 option and read out as an observational array.

Note: It would be nice to eventually allow for blocks rather than long linear transects to be sampled, as square blocks might more realistically correspond to the kind of sampling that would be done by a real observer. I don't think that this would make too much difference in terms of finding sampling error, as there is no bias to resources movement in one direction; hence, the turnover of resources for any particular number of cells will be the same for any N cells sampled. It also stands to reason that this error should be normally distributed as the number of sampling attempts becomes large, and the error should be mean centred around the actual population size, since the probability of missing and double counting would seem to cancel out exactly. This might eventually lead to analytical estimate of observation and error actually being reasonable under some conditions.

It will be helpful to link the appropriate element of the interaction array (Jacobian matrix) to the actions in the landscape_actions function in user.c. As of now, the amount of increase in crop yield (and decrease) is hard-coded in the function, but it really should be linked with the appropriate diagonal element in the interaction array -- increasing or decreasing a cell's value by the magnitude in the array element.

Apparently found a strange plotting issue given a specific set of parameters, where 'actions made' plot seems to plot <0, whereas action array definitely does not show this.

Reproducable example:

library(GMSE)
sim <- gmse(land_ownership = FALSE, stakeholders = 50, observe_type = 2, agent_view = 0.5,
            res_death_K = 2000, manage_target = 1200, RESOURCE_ini = 1200, 
            user_budget = 100, manager_budget = 1000, res_consume = 0, 
            plotting = TRUE, time_max=10)

# Note negative values plotted in bottom right panel:
plot_gmse_results(res = sim$resource, obs = sim$observation, land = sim$land, agents = sim$agents, paras = sim$paras, ACTION = sim$action, COST =sim$cost)

# Summarise 'action' array manually, note values not zero:
test <- lapply(sim$action, function(x) sum(x[1,9,2:50]) )
test <- unlist(test)
plot(test, type='l')

In light of the above reasoning, I think I'll plan to switch AGENTS to an int type, then see how this affects things:

The method sampling for case 0 is too confusing. Sometimes it means randomly sampled fix_mark individuals from the population, and sometimes it means sample within a particular range of view. Change this so that the switch functions have four clear cases:

0. Sampling with a range of view (i.e., don't rely on the fix_mark > 0 for switching methods)
1. Linear transect
2. Square transect
3. Sampling fix_mark times randomly on the landscape.

This will avoid a lot of hassle, even if the code for cases 0 and 3 end up looking the same, or very similar.

Need to change the procedure so that individuals don't get marked a second time if they're being seen the first time. In other words, we need to mark in the first iteration, then recapture if already marked. Might use a new column for this.

In other words, use the field_work function, checking what iteration it is, then make a new column only for initial captures (i.e., increment if itr = 1).

When edge_effect = 0, and therefore nothing happens when resources and agents move off of the edge of the landscape, R crashes. This is almost certainly due to some sort of memory leak. This is a low priority issue at the moment because I cannot think of a reason why anyone explicitly want the model to just ignore resources moving off of the landscape if someone wants something other than a torus (edge_effect = 1), such as a reflective edge or emigration upon leaving the landscape, this should be explicitly coded into the edge_effect function in utilities.c. Until someone asks for it, I'll stick with a torus.

Consider adjusting sensitivity each time step by re-inserting culls/cull_attempts into a new time step, and likewise for other variables.

It would be useful to allow for simulations to dynamically adjust the caution that the manager has when changing actions. At the moment, managers always assume that some specified number of actions will be performed by users, and this number does not change over the course of the simulation. But managers might be able to use the history of user actions to learn to be more or less cautious when setting new policy.

Need to add a case to the mover function (case 3). Currently, we have:

Case 0: no movement
Case 1: uniform movement in any x and y direction
Case 2: Poisson movement in the x then y direction

What is missing is a case in which movement happens with some Poisson probability, then that movement goes in any uniform direction (note: case 2 is a bit unrealistic -- perhaps enough so that it's worth getting rid of entirely, but save for now).

I can't tell if I'm just overloading R by running the simulation too many times too quickly (clicking to fast), or if there's actually a bug here. But when I comment out the below lines of code in the send_agents_home function of user.c, things seem fine.

while(agent_ID != landowner){
    do{
        agent_xloc = (int) floor( runif(0, xdim) ); 
    }while(agent_xloc == xdim);
    do{
        agent_yloc = (int) floor( runif(0, ydim) ); 
    }while(agent_yloc == ydim);
    landowner = (int) landscape[agent_xloc][agent_yloc][layer];
}

When I re-run the code quickly in succession, the above (I think) will very rarely crash the G-MSE program. I can't figure out why yet.

Allow agents to move in each time step, permanently, in some way. This might be best done through the anecdotal function. As of now, they go back to their original place at the end of each time step, and it would be good to have an option to let them move all around the landscape.

There are at least two typos in the most recent SI section (gmse_apply supplement). This should be fixed.

For some reason, the initial seed of the genetic algorithm appears to be having an effect that it shouldn't. When there are no individuals seeded in the genetic algorithm from the previous generation, the agents appear to go under-budget. It's not clear why this is the case. Oddly, managers appear to use a budget of 250 despite it being set at 300 given any seed greater than zero. When the seed is zero, the budget for setting costs drops to ca 100 for reasons that are not at all clear to me. For stake-holders, the cost drops to a fraction of its set budget (about a 30th of it). Yet, the stake-holder cost is still too low even when a seed of 20 is set; most stake-holders spend ca 1/6 of their budget when they should be forced to spend all of it.

For some reason, some initial testing seemed to suggest that resource population growth increases with the number of stake-holders, even if those stake-holders are hostile to the resource. Some further testing confirmed that stake-holders don't engage in actions there are more than two of them -- it's possible that I hard-coded something during testing, but it needs to be fixed. For now, I'm shifting the default testing options to 3 stake-holders to isolate the issue.

In functions in the genetic algorithm res_to_counts and policy_to_counts, the projected consequences of actions needs to be fine tuned. As of now, it predicts one fewer resource from movem, killem, and castem, and one more resource from feedem and helpem in res_to_counts. In policy_to_counts, it predicts one fewer resource for killem and one more resource for feedem and helpem. Really, there should probably at least be an option to use more precise estimates of what will happen. For the user function, this matters a bit less because stake-holders typically just want more or less of a resource. Managers, however, are trying to hit a middle ground a lot of the time; it is also more reasonable to assume that they have demographic information on the resources of interest.

The unresolved problem

As discovered by @jeremycusack -- if the function gmse_apply is called repeatedly with the same old_list, then R will crash. The code below replicates the crash.

sim1 <- gmse_apply(get_res = "Full");
sim2 <- gmse_apply(get_res = "Full", old_list = sim1);
sim3 <- gmse_apply(get_res = "Full", old_list = sim1);

The problem is a bit perplexing because the crash happens even if sim2 is removed and gc() called before running the third line above. The problem also does not occur if alternate functions are run in gmse_apply, which leads me to believe that the issue is with C. Specifically, it appears that the cause of the crash is related to the resource function.

The reason that the crash occurs is that the old_list is somehow modified in the process of running gmse_apply. Specifically, the sim$PARAS[33] is updated to the new resource model. I will figure out how to stop this from happening.

A temporary workaround fix

Here is a way to avoid the crash; it should be emphasised that this is really only needed if gmse_apply needs to be called with the exact same old_list. The crash will not occur if the old list differs among calls to gmse_apply, as typically will occur within a loop where the output is constantly being updated and the new gmse_apply call pulling results in from the most recent call. The workaround requires copying sim1 into two separate objects, saving each of those two objects, deleting all three objects, then reloading the two copied objects. Below is a successful example.

sim1  <- gmse_apply(get_res = "Full");
sim1a <- sim1;
sim1b <- sim1;
save(sim1a, sim1b, file = "sims.RData");
rm(sim1);
rm(sim1a);
rm(sim1b);
invisible(gc()); # Garbage collect to remove all junk
load("sims.RData");
sim2 <- gmse_apply(get_res = "Full", old_list = sim1a);
sim3 <- gmse_apply(get_res = "Full", old_list = sim1b);

Now we have successfully created a sim2 and sim3 run from the same old_list. Note that we can recover sim1 by just defining sim1 <- sim1a.

Of course, it will be easy to make this kind of transect sampling random instead of comprehensive over the landscape. This can be done by simply randomly choosing the positions of block on a landscape some obs_iter of times. This could allow an estimate of population size by considering density (i.e., assume that the number counted in a sampled block reflects the density of the larger landscape of known size), as was done by Nuno et al. (2013). This shouldn't take much time to code and test.

As noticed by @jeremycusack -- which we both thought was a bug -- the gmse_apply function requires re-stating custom functions as arguments even when such arguments are available in the old_list argument. It shouldn't be difficult to get gmse_apply to recognise when it has a custom argument from an old list, then use that argument (overriding the default that is assumed when the argument is not specified).

Resources are retaining their values of helpem and feedem after being helped for one generation. Worse, they are passing their inherited characteristics on to their offspring. This needs to be changed so that agent actions have the temporary effect of increasing offspring survival probability or reproduction -- else populations will never run the risk of crashing.

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 (m). This is simply calculated below:

(1 - m)^A = 0.5

The above can be re-arranged to find A,

A = log(0.5) / log(1 - m).

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

The value of m 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 m is simply equal to remov_pr. If res_death_type = 2 (density dependent), then m could be found perhaps using something like the following:

m = (RESOURCE_ini * lambda) / (RESOURCE_ini + RESOURCE_ini * lambda)

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).

The resource-wide parameter values (e.g., carrying capacities, movement types) will need to be either:

1. Defined at the individual scale so that each individual resource has it's own value which can then be called in the resource function as necessary, and/or
2. Input as vector in the base gmse function, the length of which could determine how many times resource is called in one time step (one for each type of resource, potentially, if carrying capacity is type specific -- or carrying capacity could be applied within a type in c -- perhaps more efficient, but would require to read in multiple K somehow, either through the paras vector or in the resources array -- or something else. How to do this best will need to consider both computational efficiency and clarity/ease of coding.

Note that:

• res_remove can already be called in a type-specific way by resource, so it might just be better to call resource once and somehow input variable numbers of K into c. I'll need to think more about this, but it could be something like assigning each individual a competition coefficient alpha for how it is affected by each other type of individual. Intra-type competition could then be modelled generally, with K defined by its inverse. Meanwhile, inter-type competition coefficients could also be useful.

• Along these lines, it's also worth considering an option allowing only one resource per cell (equating to a local alpha and K of one). This might be worth making its own issue later.

• If we were to call resource multiple times, we would also need to paste arrays together in R. This wouldn't be terrible, but it could lose some efficiency unnecessarily, and I don't see the benefit.

When only one stakeholder is included gmse_apply, R crashes. Include a simple error message when this is set.

The case 0 observation type is consistently underestimating the true population size. This could be caused by a calculation that assumes that the size of the sampled area is larger than it actually is, or that the size of the landscape is smaller than it actually is; either way, the observation.c file needs to be double-checked and potentially debugged.

Currently, the predicted effects of a manager's actions are set to values that, heuristically, appear to work in the genetic algorithm. This is adjusted with the manager_sense parameter, which has a default of 0.1, such that the manager assumes that if they set costs to increase culling by 100 percent, it will actually only increase by 10 percent (as not all users are going to necessarily cull if given the opportunity). Like real-world management, this is heuristic and results in uncertainty, but future versions of GMSE could dynamically modify this value during the course of the simulation based on real knowledge of how policy changes have affected user actions in previous time steps.

Currently, more than one resource type is permitted, but this is not offered/visible to users of the software. A next major version of GMSE could have multiple resource types with resources actually interacting with one another (could borrow future development code from EcoEdu). Simple interactions could include competition and predator-prey functions in the resource model. The code is also already ready for managers and users to consider multiple resources in making policy decisions and actions, respectively.

The action array has three columns of util, u_loc, and u_land, which represent the utility of a resource, whether or not actions on the resource are restricted to the user's land, and whether or not the utility of the resource is dependent on it being on the user's land. Currently, any positive values correspond to some cost in the cost array, which means that they are changed to zero when the cost is high. In essence, these three columns represent identity, while the remaining rows to the right represent actions. Ideally, we don't want the users to be affecting, or the constrain_costs function changing, util, u_loc, and u_land columns -- only the ones to the right.

What needs to happen next is for util, u_loc, and u_land columns to be untouchable by the genetic algorithm when the first row in the action array (agent) is negative -- corresponding to direct actions of the user on resources or landscape layers. Remaining util, u_loc, and u_land should be touchable. Hence, within constrain_costs, it is necessary to block adjustment to the relevant columns.

Because resource() might need to be called multiple times for each type of agent, it would be useful to allow for calling the C function multiple times in R. To do this, it might be useful to have the type input as a variable or a vector in the resource() function in resource.R, and use an if statement to handle differently if a vector.

If a vector, a NULL could be created, then a for loop could cycle through all agent types, iteratively augmenting a data frame. This might not be ideal, as the for loop could add more time than desired if it's done in R.

Perhaps another idea would be to allow (or force) the type to be read into c as a vector, with the loop going through the resource.c file -- the reason for not doing this was to allow different types of observation, but maybe it would be worth it to allow all the relevant variables be read in as vectors? That way everything could be done in c. Still, if, e.g., each individual were observing separately, this would be a lot to read into c.

Maybe a compromise between the two is better -- by allowing multiple agent types to be passed to c, but not multiple observation types, the observation function would only need to be called a few times, at most -- maybe for the manager, and for each type of stake-holder in the same way.

Issues with the below:

==45075== at 0x167F488A: res_to_counts (packages/tests-vg/GMSE/src/game.c:265)

and

==45075== at 0x167FD2A5: user (packages/tests-vg/GMSE/src/user.c:654)

The cast should be explicit in each case -- though I'm not sure why the first needs a cast to (double). Resolve, then re-run valgrind to double-check.

Originally, I had the idea to use the global vector para as a way of storing information easily and using it across all of the models. The vector para would store key information about pretty much everything, then be dynamically updated as need be from higher level functions in the model. In the last two months of coding, I have been specifying parameter names in functions explicitly, which has made sense during the coding process for my own writing, but it will be beneficial to clean all of this up later by reading para into these sub-functions that otherwise have sometimes about a dozen arguments. Most functions would then have considerably fewer arguments, and the description of variables stored as vector elements in para could be immediately defined within sub-functions and used by name thereafter. The whole program would then have a similar feel of reading in key arrays and vectors and then specifying the key variables within sub-functions.

Need to test the code for very large population sizes.

Currently, estimating total population size using a sub-sample of observed area and assuming that the density of this sub-sample reflects global density (method = case 0) only works when one sub-sample is taken. There are multiple ways of fixing this so that the population size estimate takes into account multiple sub-samples. It would be a good idea to think about the most efficient way to do this and program it into R (perhaps with tapply to start, but eventually in the manager.c function, maybe).

It would be useful to incorporate observation error into the simulations more directly. This could be affected by one or more variables attached to each agent, which would potentially cause the mis-identification (e.g., incorrect return of seeme) or mis-labelling (incorrect traits read into the observation array) of resources. This could be done in either of two ways:

1. Cause the errors to happen in 'real time' -- that is, while the observations are happening in the simulation. This would probably be slightly inefficient, but have the benefit of being able to assign errors specifically to agents more directly.

2. Wait until the resource_array is marked in the observation function, then introduce errors to the array itself, including errors to whether or not resources are recorded and what their trait values are. These errors would then be read into the obs_array, which is returned by the function.

Add multiple dimensions to the landscape, allowing different layers of landscape properties -- probably as a list in R, but a 3D array in C (or multiple 2D arrays?)

Currently, the observation model only records resources into the observation array if they are of a particular type1, which is specified in the para vector and used to produce a data array with only one type of specified resource. Originally, this seemed like a good idea, but after spending some initial time writing the management model, I don't think there is any need nor good reason to restrict observation to a specific resource type. Instead, all types should be marked and moved to the observation array. Then, if management analysis wants statistics for only one type of resource, its very easy to use an if statement to check that the type is appropriate. It's much easier to ignore parts of the array than to make more than one array when needed through multiple calls of the observation function.

To fix this, it shouldn't be much more than a simple removal of specifying res_type values in observation.c. When there is only one resource type, all calculations should proceed normally, but when more resources are introduced, an if is needed for both management and plotting (different groups of resources could even be made, ignoring subdivisions, by skipping the if if the type specified to look at equals -1.

gmse_apply can crash if the manager_budget is too high.

Make sure that the estimate of the population as reflected by manage_freq is accurate in gmse_table.

The gmse argument converge_crit applies a convergence criteria to the genetic algorithm such that mean fitness increase needs to be below some magnitude to terminate one genetic algorithm run. At the moment, this only applies to users in practice because the manager's fitness is assessed by minimising its deviation such that each individual in the POPULATION has its deviation from utility subtracted from the maximum individual's deviation.

The number of landscape layers in the COST and ACTION arrays is too many -- the utility_layer adds one in for every unique resource instead of for every unique landscape layer. This is an easy fix by adding an option to the function to specify landscape layers.

The notebook has started to have a hard time loading, potentially due to some R code inserted on 22 SEP 2017. If this persists, it's probably best to just remove (or comment out) this R code.

In the last couple versions of GMSE, I have been unable to increase the size of the paras vector due to a heap buffer overflow that CRAN finds.

ERROR: AddressSanitizer: heap-buffer-overflow on address 0x61800bc11bf8 at pc 0x7f9aa0c65caa bp 0x7ffcfc507a90 sp 0x7ffcfc507a88

The heap-buffer-overflow referred to the following line of code:

0x7f9aa0c65ca9 in new_act /data/gannet/ripley/R/packages/tests-clang-SAN/GMSE/src/game.c:466:2

This is going to be a restriction in future development, so I am flagging this up as an issue to address. I think the issue is this line (equivalent exists in the other submodels). I should just initialise a pointer for paras (e.g., paras_ptr), then allocate memory to paras and read in all of the values. This should make things consistent and avoid any surprises in future CRAN updates.

A suggested enhancement by @jeremycusack is to get manage_freq to work in gmse_apply and effectively skip the management step when old_list is not NULL. The way that I would do this is to have an if statement in the main gmse_apply function that checks for the input manage_freq = -1 (negative one being easy to interpret). If that is the case and there is an old list with user costs in place, then management could be skipped, allowing for the manage_freq option currently restricted to the gmse function to also be used in gmse_apply.