bradduthie / resevol Goto Github PK

Valgrind returns a strange segfault traced to the parents.c file, but only when reproduction is asexual and there is more than one pesticide. Example below.

library(resevol);
set.seed(Sys.time());
gmt       <- matrix(data = 0, nrow = 2, ncol = 2);
diag(gmt) <- 1;
mg        <- mine_gmatrix(gmatrix = gmt, loci = 4, layers = 2, indivs = 100, 
                          npsize = 100, max_gen = 2, prnt_out = FALSE);
sim       <- run_farm_sim(mine_output = mg, N = 100, xdim = 40, ydim = 40, 
                          repro = "asexual", time_steps = 3, pesticide_rotation_type = 1,
                          print_inds = FALSE, print_gens = FALSE, pesticide_number = 2,
                          print_last = FALSE, get_stats = FALSE);

This produces the below with valgrind.

==12609== Invalid read of size 8
==12609==    at 0xF933AEF: add_asexual (parents.c:785)
==12609==    by 0xF933EC1: make_offspring (parents.c:897)
==12609==    by 0xF9350EA: sim_farming (sim_farming.c:166)
==12609==    by 0x4F380CF: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F38605: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F792FE: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F835DF: Rf_eval (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F853FE: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F861C6: Rf_applyClosure (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F7A26B: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F835DF: Rf_eval (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F853FE: ??? (in /usr/lib/R/lib/libR.so)
==12609==  Address 0xacd7748 is 0 bytes after a block of size 344 alloc'd
==12609==    at 0x4C31B0F: malloc (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==12609==    by 0xF93508F: sim_farming (sim_farming.c:160)
==12609==    by 0x4F380CF: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F38605: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F792FE: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F835DF: Rf_eval (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F853FE: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F861C6: Rf_applyClosure (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F7A26B: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F835DF: Rf_eval (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F853FE: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F861C6: Rf_applyClosure (in /usr/lib/R/lib/libR.so)
==12609== 
==12609== Invalid write of size 8
==12609==    at 0xF933B0E: add_asexual (parents.c:788)
==12609==    by 0xF933EC1: make_offspring (parents.c:897)
==12609==    by 0xF9350EA: sim_farming (sim_farming.c:166)
==12609==    by 0x4F380CF: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F38605: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F792FE: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F835DF: Rf_eval (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F853FE: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F861C6: Rf_applyClosure (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F7A26B: ??? (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F835DF: Rf_eval (in /usr/lib/R/lib/libR.so)
==12609==    by 0x4F853FE: ??? (in /usr/lib/R/lib/libR.so)
==12609==  Address 0x0 is not stack'd, malloc'd or (recently) free'd
==12609== 

This makes no sense. I can't figure out why changing pesticide number should affect the reproduction of only asexual individuals.

• Crop and Pathogen info?
• Summary details for patch? (diversity, ressitance and pop. size)

http://r.rosemckeon.co.uk/helicoverpa/

• Crops
• Pathogens

Sidebar should toggle in and out of view, and patches should resize to fill all the remaining window space properly.

If this isn't done, then the oldest will feed first, causing demographic effects that don't make sense.

After looking through some simulation results, I realised that the landscape is not being cleared, so pesticide just goes down everywhere and never comes off with the landscape change. This is easy enough to do by creating separate clean_pesticide and clean_crop functions to replace landscape values with zeros each time the landscape changes. Then both the crop and the landscape should be removed and replenished when the land changes for each respectively.

I'm not sure why individuals were able to keep feeding before the landscape changed, as the crop amount should be going down. Need to check on this.

In the next version of the software, it would be good to allow users to customise file names for the different CSV outputs. Rather than try to read the names into C, I think completing a simulation with an if statement (for non-default) and file.rename to change the default would be the best way to do it.

Loop everything through a while loop by getting the living individuals back into a pest array at the end of the time step before starting the next.

• Diversity
• Population size
• Resistance

Threshold values could influence the style of the summary - ie: text could be green when the values are good, or we could have tick icons etc.

There really should be column headers in the last time step individual-level CSV files. The columns can be worked out by looking at the code, but I don't think that people should be expected to do this. The tricky bit is that there will need to be some kind of clever way to assign the headers based on the model parameter values, but this shouldn't be too difficult. The first 100 or so columns are the same, and the rest just depend on the genome structure. This can be done in R first before any of the printing in C.

I really like the intuitive graphic and don't want to mess things up. So subject to avoiding breaking any of the messaging, it would be cool if there were a few extra features to illustrate the outcome of the modelling. Some suggestions, with no particular order, and with no particular urgency:

1. Instead of just illustrating 100 larvae, sample a proportion of the survivors so that the count more accurately reflects density.
2. Can we separately illustrate fitness in the two dimensions of the landscape, e.g., with size indicating crop eating ability, and maybe x for eyes if the caterpillar dies from infection (or something better).
3. Could the fields be coloured and/or textured according to the crop and pathogen combination that was last grown on them? This might require a legend to interpret the backgrounds, but would be useful for illustrating differences across the landscape?

No worries if you disagree or any of this is not sensible, of course....

A simple function is needed so that individuals increase their age within each time step.

There seems to be a memory issue in v0.2.0.8 that only applies to Windows systems (CRAN check). I'm not quite sure what the issue is here, but it will be necessary to use a Windows system to isolate the problem and resolve. WinBuilder doesn't pick it up.

Hi All, From some constructive feedback by @MattTinsleyStirling -- I am considering adding a supporting information (SI) document (and package vignette) that provides suggestions for how to model a type of system or situation using the resevol package. That is, if someone said, "How would I model pesticide resistance to GM crops in an annual system?" (or something slightly more specific), the SI would make a suggestion for how to do it. What would be massively helpful here would be if you could:

1. Comment on whether or not you think that this is important for the initial submission to Methods in Ecology and Evolution, or if this would be something to add later to the package.
2. Suggest some systems to model. E.g., can you reply to this asking "How would I model X", where X is a specific system or scenario that you think might be useful to model.

I do like the idea of having a document that illustrates how these abstract options in the software could usefully be applied to modelling specific systems. It might minimise the risk that readers miss an opportunity to model their system because it's not clear enough how a particular set of options could be applied.

We can use this issue to propose and discuss potential ideas for simulations to address interesting scientific questions.

I think it would be good to include the option to input some kind of vector or matrix as an argument in run_farm_sim that would specify the transition from one crop (or pesticide) to the next. This could be something as simple as a vector that specifies the order of each crop or pesticide. The challenge would be to implement some pattern where one crop repeats before another (e.g., 1, 2, 2, 3, 1, 2, 2, 3). A cell could perhaps cycle through the vector, with different initial starting points among cells?

The pests need to accumulate biopesticide when they land on a cell, assuming the correct age. The details of consuming the biopesticide can be written almost identically to the feed function, because the general idea is the same. We need to add a column in the pest array for summed biopesticide consumption, but this is fairly easy.

I still need to write reproduction that is consumption based.

Need to rotate the landscape values, or have the option to do so, given something like ts %% X == 0 where ts is the time step and X is the frequency of landscape rotation. The function to do this needs some options, the easiest of which will be to randomly reassign something for each unique block.

I'm going to add in sexual reproduction, which makes things challenging if we need to have diploid organisms, but should work with the current setup. Because the network is fixed from loci to traits to have a predefined trait covariance given standard random normals, the variation of the allele values needs to be preserved with diploidy. This is simple enough to do by initialising the allele values with a variance of 1/sqrt(2), then making the allele effects additive (i.e., a1 from mother plus a2 from father equals the offspring allele value). We don't need to control the variance for the whole loci to offspring network in this way, even if we want that network to mutate, because its values are selected through the evolutionary algorithm. We just need to either divide them by 2 so that each copy sums to the desired value (probably easiest), or take the mean when there is evolution of the G matrix. I think summing is computationally most efficient and consistent with the loci.

In sum, all individuals will get two sets of alleles (diploid) with additive effects, initialised from rnorm(sd = 1/sqrt(2)), and two sets of network values initialised by val * 0.5.

Note, this adds an extra step where the values need to be summed before producing the matrices multiplied to get the traits, but this might not take too much extra time because the values need to be placed in the arrays through a loop anyway for each newly initialised individual.

It appears that some individuals are being included in the summary statistics if they are dead but not yet removed from the population. This isn't necessarily a bad thing (they'll be removed in the subsequent time step), but maybe best to just not include these individuals when calculating summary statistics for a time step.

There needs to be a function for death according to some different death types. The most simple way to do this is to cause death if not enough total food is consumed or too much biopesticide is consumed. This might make separate columns for needed food and maximum biopesticide of different types. Note that there needs to be an age set for each too, but perhaps this can be dealt with later. One issue is whether or not to wipe food accumulated each time step -- seems like there should be an option for this (easy to add to paras).

So far interface patches are built with made up output in the following format...

{
"crops": [1,2],
"pathogens": [0,1],
"helicoverpa": [
{"crop":1, "pathogen": 0, "genotype":1, "resistance": 30},
{"crop":2, "pathogen": 0, "genotype":1, "resistance": 30}
]
}

This is an array of the crop species as numbers for each patch, and the same for pathogens. The size of these arrays defines the number and arrangement of patches.

Then the Helicoverpa array contains every individual. I've used the crop and pathogen values for each individual to place them in the correct patch -will that always work? Maybe we just need a horizontal and vertical patch id value instead so we can give each individual a coordinate like placement...

{
"crops": [1,2],
"crop_ids": [1,2] ,
"pathogens": [0,1],
"pathogen_ids" : [1,2],
"helicoverpa": [
{"c_id": 1, "p_id":1, "genotype":1, "resistance": 0.8}
] 
}

The individuals also have genotype and resistance values to define their colour and size. The types of numbers you give me can be whatever size you think best - I can fiddle with the styling based on those numbers after.

From Brad:

Just looking through the current output now and trying to get a handle on the json. Sorry if I'm missing the point here somewhere! Would the following output format work?

{
    "crops": [1, 2, 3],
    "pathogens": [1, 2, 3],
    "helicoverpa": [
        {"crop":1, "pathogen": 1, "xloc":1, "yloc":2, "path_allele_1":1, "path_allele_2":2, "crop_allele_1":1, "crop_allele_2":1},
        {"crop":2, "pathogen": 1, "xloc":2, "yloc":2, "path_allele_1":2, "path_allele_2":2, "crop_allele_1":3, "crop_allele_2":1}
        ]
}

I assume that the top "crops" and "pathogens" is the set of crops and pathogens available on the landscape? I'm totally happy to fill this in any way that works best for you -- here's what information we have:

LAND: pathogen and crop on each cell in location x and y.
PEST: Each individual pest: ID, sex, x-loc, y-loc, 2 pathogen alleles, 2 crop alleles.

Sorry, I think I'm just confused as to how the information produced in the toy matches with the json format!

I don't know if this is a problem or not, but in writing up the documentation, I recognised that the food consumption that directly translates to survival and reproduction is never 'cleared' during an individuals lifetime. In other words, if an individual consumed 4 units of food and the threshold for surviving was 1 and the reproducing was 2, then each subsequent time step the individual would survive and produce 4/2 offspring. That is, the consumption increased per offspring production each time step rather than once, so an individual that continued to survive for many time steps would get to also produce more offspring consistently in each time step. This isn't necessarily bad, but I want the option to change it. I think that the change will be useful @rosemckeon @luc-bussiere @RosieMangan @MattTinsleyStirling

Worse, the food accrued from previous time steps would continue to build up -- as in, if an individual ate 1 unit of food two time steps ago, then they would get credit for eating 1 unit again on top of whatever they eat in the current time step. Obviously this is just wrong. That has been cleared with a new function refresh_consumed in the feeding.c file. All tested and appears to work.

Next, I introduce a metabolism into the model, which decreases the food consumed by an individual by some real value in each time step. The value can also be a trait, so some individuals can burn food faster than others (which could potentially be interesting to vary). This will be committed and pushed soon, but it has been tested and appears to be working okay. Annoying that I didn't include this sooner, but glad that I found it when writing the documentation.

I'll just close the issue once everything is pushed and double-checked.

Diploid simulations are taking too long; there is no reason that I can see why they should be taking orders of magnitude longer than haploid simulations. I need to figure out why this is and fix it as best as possible. There's just no reason that I can tell why add_sexual, sire_genes, or assign_sire should be taking so long. The way that the traits are inserted, mutation occurs, or inbreeding coefficients are calculated also just shouldn't take that much time. Maybe calling is_in_range so many times is the problem?

It would be useful to write a function that takes the output of mine_gmatrix and explores as wide a range of loci values as possible (e.g., four standard deviations around 0 for each locus) in all possible combinations. This would seem to give the possible constraints of trait space rather than the constraints inherent to the initial genome.

Patches should use all available screen space and Masonry should arrange individual Helicoverpa. FontAwesome can be used to control the happiness of individuals based on resistance.

For some reason, after mining for a gmatrix, certain mines cause the simulation take much longer to run than others. Some mines cause simulations of 100+ generations to be run in under a second, whereas others take nearly a minute for some reason. It is worth printing off the timestamps for each function in the sim_farming C function to see what the rate limiting step is, then trying to match that with the mine gmatrix details.

After a recent lab meeting, I think it would be interesting to introduce haplodiploidy as an option in the simulations. We could then use the package to answer some interesting questions about inbreeding in diploid versus haplodiploid systems.

For both sexual and asexual reproduction, in the count_offspring function, there needs to be some if statement that only allows individuals to reproduce if they are within an acceptable age range as determined by their individual column values. The mate_in_range function also needs to include an if statement to avoid having individuals identify a mate that is not within an acceptable range.

For sexual reproduction, the sire assigned in the assign_sire function also needs to include an if statement as in the mate_in_range function so that the sire is only assigned if an acceptable age.

New code is needed to run multiple iterations of the toy_simulate_resistance function across different parameter values. The key issue here is to get a range of parameter space for the heterogeneity of the landscape, which I plan to do by (at first) running the same set of parameter values, but building a new function run within the toy_initialise_land function. Currently, crops and pathogens are assigned to cells randomly, but some sort of toy_set_land_block could ensure that the whole landscape is separated into $N$ discrete blocks. This line of code from GMSE might be helpful to set blocks.

Once the blocks are set, we can then change $N$ to increase or decrease the scale at which heterogeneity of crops and fungal biopesticides are applied (i.e., high $N$ means more continuous blocks are created, low $N$ means few blocks are created -- N < 2 means that there's only one big block). Crop and fungal biopesticides would then be assigned based on some simple rules.

As suggested by @MattTinsleyStirling -- write a function to run a sensitivity analysis for the effect of each locus on each trait. This should be a fairly straightforward task of changing the value of each locus by an increment of some kind (or maybe multiple increments at different intervals in each direction), then getting the partial change in traits as a response. This could be an intermediate function between mining the gmatrix and simulating, e.g., gmatrix_sensitivity.

A function that moves individuals, ideally in a random order, needs to be written. There should also be an option that allows individuals to eat each time they move, and perhaps settle given acceptable food and mate given a partner on the same location.

• Feeding should first be allowed for all living individuals of an acceptable age.
• Movement should then follow in a new function: given more 'bouts' of movement, an individual will have the opportunity to land on a cell and feed before moving to a new cell (all in a random order). This will allow simulations where individuals move around and feed if desired, while also restricting simulations to a simpler feed then move once scenario by setting the 'bouts' parameter to 1.
• It would be nice if all individuals could move around and females could mate upon landing on a cell with a male nearby, but this would require some kind of record of all the males with which each female crossed paths. This is certainly possible, e.g., with an N by N binary array that indicates overlap in a time step, but this would then require a rewrite a new option in the assign_sire function to perform the lookup. This sounds like a lot of fun, but maybe something to put off for later.
• Note that the order of mating can be controlled by defining the age window in which mating and movement are allowed. Time steps are not quite to be thought of as seasons here.
• Death needs to be modelled as a function of not getting enough to eat or getting too much pesticide. Independent mortality should also be an option (or perhaps simple density dependent).

The next step is to sort out feeding and movement though.

A good point raised by @MattTinsleyStirling is that there is no option for a cells crop value to increase over time following crop rotation. This made me think that it might be easy to add such an option . The default value could be zero, but nonzero values could specify the proportion to add in each time step. For example, if there are two crops and we set the argument crop_growth = c(0.25, 0.5), then in each time step, the cell values for crop 1 could increase by 25 per cent, and the cell values for crop 2 could increase by 50 per cent. I actually don' t think that this would take too long to code (few hours, maybe?). Is it worth doing, you think, before submitting?