This notebook tracks progress on the development of resevol R package, a BBSRC FAPESP Newton Funded Project to model how the spatial scale of heterogeneity in fungal isolate application and crop plant cultivation affects variation in the selection landscape for biopesticide resistance.


Go to comments

Contents:

Project updates

Project comments


Project updates:

Update: 30 JUL 2022

I don’t know why I didn’t think of this sooner, but there really should be a feature that allows users to set the initial mean values of traits. This is easy because trait means are of course independent of the covariance structure of the traits, so there just needs to be some specified increment to be added to each trait throughout the simulation. I have introduced this in a new version 0.3.2.0 with the trait_means argument, which simply takes a vector of the means for each customised trait. By default, values are set to zero (i.e., what they always were originally), but the argument can also take a vector of the same length as the number of traits, with each element position corresponding to the specific trait. That is, trait_means[1] is the mean value of "T1", trait_means[2] is the mean value of "T2", et cetera. An error is given if the user attempts to specify a trait_means vector that has a length different from the number of traits. Just to show this in action, I will use the resevol mine_gmatrix output mg_v1, which has the covariance structure below.

##            [,1]       [,2]       [,3]       [,4]
## [1,]  0.9936103 -0.5106387  0.2324040  0.2367691
## [2,] -0.5106387  1.0193888  0.2573216  0.2723028
## [3,]  0.2324040  0.2573216  0.9334559 -0.4480018
## [4,]  0.2367691  0.2723028 -0.4480018  0.9689539

We have four traits, which I will use to define consumption rate of crops 1 and 2 (T1 and T2, respectively), and the consumption rates of pesticides 1 and 2 (T3 and T4, respectively). I will set the trait means to T1 = 1, T2 = 2, T3 = 0.2, and T4 = 0.2 and run a simulation with the code below.

sim       <- run_farm_sim(mine_output = mg_v1, N = 1000, xdim = 10, ydim = 10, 
                          repro = "asexual", time_steps = 20, farms = 12,
                          crop_number             = 2,
                          crop_init               = "random",
                          crop_rotation_type      = 2,
                          pesticide_number        = 2,
                          pesticide_init          = "random",
                          pesticide_rotation_type = 2,
                          food_consume            = c("T1", "T2"),
                          pesticide_consume       = c("T3", "T4"),
                          trait_means             = c(1, 2, 0.2, 0.2),
                          food_needed_surv        = 1,
                          food_needed_repr        = 3,
                          print_last              = TRUE);

At the end of the simulation (after some evolution has occurred), we can look at the evolution of traits and final trait distributions to confirm that they are around what we expect. Here is the evolution of the different traits over time.

Here are the distributions of each trait, individually.

Finally, here is a pairs plot showing the trait correlations.

Overall, this is doing what we want it to. I have also written a new test into testthat for the package. New documentation for the trait_means argument has been written, and I will incorporate this option into the working manuscript.

Update: 19 JUL 2022

Version 0.3.0.1 has been submitted to CRAN and is passing all checks. I am now a bit under halfway through writing a first draft of the software note. This note be submitted as a pre-print, and for publication, and will explain how to use the package. It will also replace the current Get Started article and CRAN vignette. The structure of the publication introduces some key background, the mine_gmatrix function and its underlying logic, the ecology and evolution of pests in the model, and a simple example. I think that two supporting information documents are also in order. The first will explain the evolutionary algorithm in more detail and basically be a version of the current vignette with maybe some minor revisions. The second will be a more advanced example with complex dynamics and traits, presenting simulations like the ones shown in this journal on 26 AUG 2021. My hope is that this plus the documentation is sufficient for anyone to start using the package.

Having submitted version 0.3.0.1, I almost immediately discovered the need for an update that is now version 0.3.1.0. I will hold off on submitting this one to CRAN just yet, as the update is quite minor. The cause was that the mine_gmatrix function kept giving a discrepancy between the stress threshold that caused it to terminate and the actual predicted stress of the network it found. The reason for this was that the algorithm was terminating due to finding unusually low stress caused sampling error for the standard random normal loci. In other words, a luck draw for random loci often caused networks to look better at producing desired trait covariances than they actually were. I’ve addressed this now by first using the mean stress for a whole population of networks to test the termination criteria (term_cri) in the mine_gmatrix function, rather than the best stress out of all networks in the population. But the mine_gmatrix function still returns the network that produces the lowest found stress.

Second, I have introduced a new diagnostic function to give the distribution of stress predicted from loci.

stress_test <- function(mine_output, indivs = 1000, reps = 10){
    
    target_VCV     <- mine_output[[2]];
    Loci_to_Traits <- mine_output[[5]];
    loci           <- mine_output[[1]][1];
    
    Ln_stress <- rep(x = NA, times = reps);
    while(reps > 0){
        sim_loci        <- rnorm(n = indivs * loci);
        individuals     <- matrix(data = sim_loci, nrow = indivs, ncol = loci);
        trait_vals      <- individuals %*% Loci_to_Traits;
        VCV             <- cov(trait_vals);
        stress          <- mean((target_VCV - VCV) * (target_VCV - VCV));
        Ln_stress[reps] <- log(stress);
        reps            <- reps - 1;
    }
    
    return(Ln_stress);
}

This simple function now forms part of the package itself, so users can run the output of mine_gmatrix and get a better overall picture on how the chosen network performs. The example used in the documentation is below (I’ve increased the replicates using the argument reps = 1000; the default is 10).

gmt       <- matrix(data = 0, nrow = 4, ncol = 4);
diag(gmt) <- 1;
mg        <- mine_gmatrix(gmatrix = gmt, loci = 4, layers = 3, indivs = 100, 
                          npsize = 100, max_gen = 2, prnt_out = FALSE);
stresses  <- stress_test(mine_output = mg, reps = 1000);

The output stresses is just a vector of stresses for some number (indivs) of initialised genomes. The histogram below shows the distribution for the above example.

This should at least make it easier for users to decide if they do or do not like the output of mine_gmatrix for use in a new set of simulations. The above example shows stress that is quite consistent, but if we let the simulation run for longer to find a better solution, the stress variance seems to increase (note that now max_gen = 400 below).

gmt       <- matrix(data = 0, nrow = 4, ncol = 4);
diag(gmt) <- 1;
mg        <- mine_gmatrix(gmatrix = gmt, loci = 4, layers = 3, indivs = 100, 
                          npsize = 100, max_gen = 400, prnt_out = FALSE);
stresses  <- stress_test(mine_output = mg, reps = 1000);
hist(stresses, xlab = "Stress values for initialised genomes", 
     main = "", breaks = 20);

It might be good to recommend this kind of diagnostic for users, particularly when simulating very complex individuals. I think that I will do this for the supporting information showing an advanced example.

Update: 28 JUN 2022

Custom crop and pesticide initialisation and rotation

I have now introduced custom crop and pesticide initialisation and rotation, as I set out to do in the 13 JUN 2022 notes below. The code appears to be stable, and I have introduced two new automated tests in the package to make sure that the feature remains stable. Below is an example of how to have a custom rotation regime for crops and pesticides; I will start with crops. First, I will use the mg_n1 gmatrix that comes with the resevol package (I could have also made a dummy one; the details are not important here).

data("gmatrices");

Now I will run a simulation. The run_farm_sim takes some simple options for crop_init and crop_rotation_type. This includes crop_init = "random", which randomly assigns each farm to one of the potential crops. There are also three crop_rotation_type arguments, including 1 (no crop rotation), 2 (random crop choice whenever it is time to change), and 3 (cycling through crops by number; i.e., 1, 2, 3, 1, 2, 3). I might add more preset options like this later, but here is how to run a simulation with two crops, random crop initialisation on farms, and random rotation between crops.

sim       <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4, 
                          repro = "asexual", time_steps = 10, farms = 4,
                          crop_number        = 2,
                          crop_init          = "random",
                          crop_rotation_type = 2,
                          print_inds = FALSE, print_gens = FALSE,
                          print_last = FALSE, get_stats = FALSE);

The new feature is that crop_init and crop_rotation_type can take custom arguments, so potentially any initialisation and cycling of crop choices is possible. With crop_init, we can specify which farm gets what crop at the start of the simulation using a fector of length farms. If, for example, we want farms 1 and 4 to have crop 1, and farms 2 and 3 to have crop 2, we could run the below.

sim       <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4, 
                          repro = "asexual", time_steps = 10, farms = 4,
                          crop_number        = 2,
                          crop_init          = c(1, 2, 2, 1),
                          crop_rotation_type = 2,
                          print_inds = FALSE, print_gens = FALSE,
                          print_last = FALSE, get_stats = FALSE);

Note that the position of farms on the landscape is shown below.

##      [,1] [,2] [,3] [,4]
## [1,]    1    1    2    2
## [2,]    1    1    2    2
## [3,]    3    3    4    4
## [4,]    3    3    4    4

Hence, this landscape would start out with farms numbered 1 and 4, diagonal from one another, applying crop 1, and farms numbered 2 and 3, also diagonal, applying crop 2. Now, suppose that we want users of crop 1 to switch to crop 2 with a probability of 0.75, and retain the use of crop 1 with a probability of 0.25. Likewise, suppose that we want users of crop 2 to switch to crop 1 with a probability of 0.75, and retain the use of crop 2 with a probability of 0.25. We can now specify this with a square matrix to input in place of the number for the argument crop_rotation_type.

change_mat <- matrix(data = c(0.25, 0.75, 0.75, 0.25), nrow = 2, ncol = 2);
##      [,1] [,2]
## [1,] 0.25 0.75
## [2,] 0.75 0.25

Note that matrix elements \(M_{i,j}\) specify the probability of changing from crop \(i\) to crop \(j\), meaning that row elements must sum to 1. This can be binary too, as in the default crop_rotation_type = 3.

##      [,1] [,2]
## [1,]    0    1
## [2,]    1    0

The matrix above cycles back and forth between crops 1 and 2. Overall, this should give a lot more flexibility in terms of specifying how crops cycle over time. We can input the change_mat directly as an argument in run_farm_sim.

sim       <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4, 
                          repro = "asexual", time_steps = 10, farms = 4,
                          crop_number        = 2,
                          crop_init          = c(1, 2, 2, 1),
                          crop_rotation_type = change_mat,
                          print_inds = FALSE, print_gens = FALSE,
                          print_last = FALSE, get_stats = FALSE);

Nevertheless, the cycling is restricted to Markov chains, meaning that it cannot include cycles of crops such as 1, 1, 2, 1, 1, 2, 1, 1, etc., where the cycling repeats complex patterns. This is mainly due to the awkwardness of trying to allow for such cycles given simulations that can have high numbers of time steps (i.e., I do not think it is worth figuring out a system that allows for these kinds of cycles as an option in the code).

Pesticide rotation options work exactly the same way as crop rotation options, so I will not spend to much time on them. Here is an example of running one.

sim       <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4, 
                          repro = "asexual", time_steps = 10, farms = 4,
                          pesticide_number        = 2,
                          pesticide_init          = c(1, 2, 2, 1),
                          pesticide_rotation_type = change_mat,
                          print_inds = FALSE, print_gens = FALSE,
                          print_last = FALSE, get_stats = FALSE);

In conclusion, this new feature will allow for a very flexible range of different crop and pesticide cycling options.

Custom landscape option

I have also added an option to customise the farming landscape, also as I set out to do in the 13 JUN 2022 notes below. The customisation is on the top layer of the landscape that differentiates independent farms. What is required is a matrix that includes a sequence of natural numbers in all elements, so if there are 4 farms, then all matrix elements must be either 1, 2, 3, or 4. Beyond this requirement, there is no restriction on where different farms are placed; the do not even need to be contiguous on the landscape. Custom terrains can be inserted using the terrain argument in run_farm_sim. Below makes a custom terrain with three farms.

custom_terrain       <- matrix(data = 0, nrow = 8, ncol = 8);
custom_terrain[,1:2] <- 1;
custom_terrain[,2:6] <- 2;
custom_terrain[,7:8] <- 3;
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    1    2    2    2    2    2    3    3
## [2,]    1    2    2    2    2    2    3    3
## [3,]    1    2    2    2    2    2    3    3
## [4,]    1    2    2    2    2    2    3    3
## [5,]    1    2    2    2    2    2    3    3
## [6,]    1    2    2    2    2    2    3    3
## [7,]    1    2    2    2    2    2    3    3
## [8,]    1    2    2    2    2    2    3    3

We can insert the terrain into the landscape with the terrain argument.

sim       <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4, 
                          repro = "asexual", time_steps = 10, farms = 4,
                          terrain = custom_terrain,
                          print_inds = FALSE, print_gens = FALSE,
                          print_last = FALSE, get_stats = FALSE);

Note that a custom terrain will override the arguments farms, xdim, and ydim, so even though the above sets farms = 4, there are actually only three farms on the landscape. I am tempted to put up a warning for this, but I am actually not sure if it is needed; my hope is that anyone customising the landscape would recognise this from the documentation.

Also note that these terrain values do not necessarily need to be farms. Through the use of a custom landscape and pesticide rotation option, these cells could represent something like diversionary feeding sites or even buildings or rivers. Consider, for example, the following matrix.

complex_terrain              <- matrix(data = 0, nrow = 12, ncol = 12)
complex_terrain[1:2, 1:3]    <- 1;
complex_terrain[3:7, 1:2]    <- 1;
complex_terrain[8:12, 1:2]   <- 2;
complex_terrain[1:2, 6:9]    <- 3;
complex_terrain[3:6, 7:9]    <- 3;
complex_terrain[1:6, 10:12]  <- 4;
complex_terrain[10:12, 7]    <- 5;
complex_terrain[7:12, 8:12]  <- 5;
complex_terrain[3:12, 3:4]   <- 6;
complex_terrain[7:9, 5]      <- 6;
complex_terrain[1:2, 4:5]    <- 7;
complex_terrain[3:6, 5:6]    <- 7;
complex_terrain[7:9, 6:7]    <- 7;
complex_terrain[10:12, 5:6]  <- 7;
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,]    1    1    1    7    7    3    3    3    3     4     4     4
##  [2,]    1    1    1    7    7    3    3    3    3     4     4     4
##  [3,]    1    1    6    6    7    7    3    3    3     4     4     4
##  [4,]    1    1    6    6    7    7    3    3    3     4     4     4
##  [5,]    1    1    6    6    7    7    3    3    3     4     4     4
##  [6,]    1    1    6    6    7    7    3    3    3     4     4     4
##  [7,]    1    1    6    6    6    7    7    5    5     5     5     5
##  [8,]    2    2    6    6    6    7    7    5    5     5     5     5
##  [9,]    2    2    6    6    6    7    7    5    5     5     5     5
## [10,]    2    2    6    6    7    7    5    5    5     5     5     5
## [11,]    2    2    6    6    7    7    5    5    5     5     5     5
## [12,]    2    2    6    6    7    7    5    5    5     5     5     5

We can visualise the terrain a bit more clearly with an image below.

We can imagine the dark maroon colour as being a river (7), and the red region on the immediate left of the river as being a nature reserve (6). To get the simulations to run treating it this way, we would only need to spray a crop and pesticide that has no effect on the pests on the river (effectively making it an empty area of landscape in terms of feeding or resistance), and a pesticide but not crop in the nature reserve. First we can initialise the crop and pesticide vectors. Assume that there are 2 legitimate crops and pesticides that the actual farms (1-5) will use.

initial_crops      <- c(1, 2, 1, 2, 1, 3, 4);
initial_pesticides <- c(1, 1, 2, 2, 1, 3, 3);

We can then clarify the matrix so that there is only a cycle between 1 and 2, but never 3 and 4 for either the pesticide or the crops. First I will do the crops.

crop_cycle <- matrix(data = 0, nrow = 4, ncol = 4);
crop_cycle[1, 2] <- 1;
crop_cycle[2, 1] <- 1;
crop_cycle[3, 3] <- 1;
crop_cycle[4, 4] <- 1;
##      [,1] [,2] [,3] [,4]
## [1,]    0    1    0    0
## [2,]    1    0    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1

Now I will do the pesticide.

pesticide_cycle <- matrix(data = 0, nrow = 3, ncol = 3);
pesticide_cycle[1, 2] <- 1;
pesticide_cycle[2, 1] <- 1;
pesticide_cycle[3, 3] <- 1;
##      [,1] [,2] [,3]
## [1,]    0    1    0
## [2,]    1    0    0
## [3,]    0    0    1

So farms (1-5) will swap between crops 1 and 2, while the reserve (6) will maintain crop 3 and the river (7) will be ‘crop’ 4. And farms will swap between pesticides 1 and 2, while the reserve and river will hold pesticide 3. Then, to make crop 4 and pesticide 3 not really exist, we just make it impossible for the pests to consume these.

sim       <- run_farm_sim(mine_output = mg_n1, N = 100, 
                          repro = "asexual", time_steps = 10,
                          crop_number = 4, pesticide_number = 3,
                          terrain                 = complex_terrain,
                          food_consume            = c(0.25, 0.25, 0.25, 0),
                          pesticide_consume       = c(0.1, 0.1, 0),
                          crop_init               = initial_crops,
                          crop_rotation_type      = crop_cycle,
                          pesticide_init          = initial_pesticides,
                          pesticide_rotation_type = pesticide_cycle,
                          print_inds = FALSE, print_gens = FALSE,
                          print_last = FALSE, get_stats = FALSE);

The output shows that the pests leave the ‘food’ on the river alone (see layer 5), while the reserve (see layer 4) gets used very much.

## , , 1
## 
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,]    1    1    1    7    7    3    3    3    3     4     4     4
##  [2,]    1    1    1    7    7    3    3    3    3     4     4     4
##  [3,]    1    1    6    6    7    7    3    3    3     4     4     4
##  [4,]    1    1    6    6    7    7    3    3    3     4     4     4
##  [5,]    1    1    6    6    7    7    3    3    3     4     4     4
##  [6,]    1    1    6    6    7    7    3    3    3     4     4     4
##  [7,]    1    1    6    6    6    7    7    5    5     5     5     5
##  [8,]    2    2    6    6    6    7    7    5    5     5     5     5
##  [9,]    2    2    6    6    6    7    7    5    5     5     5     5
## [10,]    2    2    6    6    7    7    5    5    5     5     5     5
## [11,]    2    2    6    6    7    7    5    5    5     5     5     5
## [12,]    2    2    6    6    7    7    5    5    5     5     5     5
## 
## , , 2
## 
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,] 0.75 0.75 0.25    0    0    0 1.00 1.00 1.00     0     0     0
##  [2,] 1.00 0.00 0.00    0    0    0 0.75 1.00 1.00     0     0     0
##  [3,] 0.75 0.25 0.00    0    0    0 0.75 1.00 1.00     0     0     0
##  [4,] 0.75 0.00 0.00    0    0    0 1.00 1.00 1.00     0     0     0
##  [5,] 0.75 0.00 0.00    0    0    0 0.25 1.00 0.75     0     0     0
##  [6,] 0.75 0.00 0.00    0    0    0 0.00 1.00 1.00     0     0     0
##  [7,] 0.50 0.25 0.00    0    0    0 0.00 1.00 1.00     1     1     1
##  [8,] 0.00 0.00 0.00    0    0    0 0.00 0.25 0.25     1     1     1
##  [9,] 0.00 0.00 0.00    0    0    0 0.00 0.50 1.00     1     1     1
## [10,] 0.00 0.00 0.00    0    0    0 0.00 0.25 1.00     1     1     1
## [11,] 0.00 0.00 0.00    0    0    0 0.75 1.00 1.00     1     1     1
## [12,] 0.00 0.00 0.00    0    0    0 0.00 0.75 1.00     1     1     1
## 
## , , 3
## 
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,] 0.00 0.00    0    0    0    0    0    0    0     1     1     1
##  [2,] 0.00 0.00    0    0    0    0    0    0    0     1     1     1
##  [3,] 0.00 0.00    0    0    0    0    0    0    0     1     1     1
##  [4,] 0.00 0.00    0    0    0    0    0    0    0     1     1     1
##  [5,] 0.00 0.00    0    0    0    0    0    0    0     1     1     1
##  [6,] 0.00 0.00    0    0    0    0    0    0    0     1     1     1
##  [7,] 0.00 0.00    0    0    0    0    0    0    0     0     0     0
##  [8,] 1.00 0.00    0    0    0    0    0    0    0     0     0     0
##  [9,] 0.25 0.00    0    0    0    0    0    0    0     0     0     0
## [10,] 0.25 0.00    0    0    0    0    0    0    0     0     0     0
## [11,] 0.75 0.25    0    0    0    0    0    0    0     0     0     0
## [12,] 0.50 1.00    0    0    0    0    0    0    0     0     0     0
## 
## , , 4
## 
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,]    0    0  0.0    0    0    0    0    0    0     0     0     0
##  [2,]    0    0  0.0    0    0    0    0    0    0     0     0     0
##  [3,]    0    0  0.0    0    0    0    0    0    0     0     0     0
##  [4,]    0    0  0.0    0    0    0    0    0    0     0     0     0
##  [5,]    0    0  0.0    0    0    0    0    0    0     0     0     0
##  [6,]    0    0  0.0    0    0    0    0    0    0     0     0     0
##  [7,]    0    0  0.0    0    0    0    0    0    0     0     0     0
##  [8,]    0    0  0.0    0    0    0    0    0    0     0     0     0
##  [9,]    0    0  0.0    0    0    0    0    0    0     0     0     0
## [10,]    0    0  0.0    0    0    0    0    0    0     0     0     0
## [11,]    0    0  0.5    0    0    0    0    0    0     0     0     0
## [12,]    0    0  0.0    0    0    0    0    0    0     0     0     0
## 
## , , 5
## 
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,]    0    0    0    1    1    0    0    0    0     0     0     0
##  [2,]    0    0    0    1    1    0    0    0    0     0     0     0
##  [3,]    0    0    0    0    1    1    0    0    0     0     0     0
##  [4,]    0    0    0    0    1    1    0    0    0     0     0     0
##  [5,]    0    0    0    0    1    1    0    0    0     0     0     0
##  [6,]    0    0    0    0    1    1    0    0    0     0     0     0
##  [7,]    0    0    0    0    0    1    1    0    0     0     0     0
##  [8,]    0    0    0    0    0    1    1    0    0     0     0     0
##  [9,]    0    0    0    0    0    1    1    0    0     0     0     0
## [10,]    0    0    0    0    1    1    0    0    0     0     0     0
## [11,]    0    0    0    0    1    1    0    0    0     0     0     0
## [12,]    0    0    0    0    1    1    0    0    0     0     0     0

Similarly, the pests will only be affected by pesticides 1 and 2, while 3 is effectively nothing. Note that this is of course a trick for modelling landscape areas without food or pesticide, but I do not see a reason to make this more explicit in the function (unless I receive requests due to confusion, but my hope is it is clear from these notes and documentation).

Update: 13 JUN 2022

As I prepare the pre-print and software note for the resevol package, I am thinking that might be worth doing just a bit more coding so that the publication can include four more features:

  1. More explicit crop and pesticide rotation, as mentioned in Issue #43. I think that this would be fairly easy to implement and worth the time. Currently, the run_farm_sim argument crop_rotation_type accepts either 1 or a 2. A 1 results in no rotation of crops on the landscape, and a 2 results in a random transition between types. I think that a custom transition matrix should be possible to insert, whereby each element should specify the probability of transitioning from crop i (row) to crop j (column). In the simplest case (e.g., could be crop_rotation_type = 3), this could just be a binary matrix that forms a cyclical graph. The current options 1 and 2 could even be subsumed in the transition matrix.
  2. Allow custom landscapes. I think that this should be possible by allowing users to put their own matrix in a function, with numbers on the matrix representing different farms (integers greater than zeros) or public land (zeros). A new custom_landscape function could check this matrix, then use it as the top layer of an array.
  3. Custom plotting functions. This one is tricky given the flexibility of simulations and therefore the complexity of the output, but it might be worth including at least some basic plotting functions, or summary functions.
  4. Simulated annealing algorithm for the mine_gmatrix option. This could be fairly easy to implement, and for some cases might actually find genome solutions faster than the evolutionary algorithm.

My hope is that these different tasks could be done fairly quickly, and that they will make writing up the paper go more smoothly.

Update: 06 JAN 2022

There had for some time been an annoying issue of resevol CRAN checks having ERRORs on Windows systems with ix86 architecture. As it turns out, this was caused by some errors in specifying the memory for arrays, which was fixed in this commit. Just to make sure that everything is clear and correct, I have made the cast much clearer in a new commit. Tests using the rhub package and on Winbuilder confirm that everything works as intended. Valgrind found no memory leaks; everything appears fine, so I have updated the version to v0.2.0.9 and submitted the fix to CRAN.

Update: 26 AUG 2021

I have now run the same analysis that I did on 19 AUG 2021, but with traits that evolve independently. This was done by setting a covariance matrix to be the identity matrix for four traits. Here is the mine_gmatrix output that was used for the simulations.

## [[1]]
##  [1]    12.00     6.00  2000.00 18000.00     0.01     0.01  3200.00     0.01
##  [9]  1200.00     6.00    -9.00     0.10     0.00
## 
## [[2]]
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1
## 
## [[3]]
##              [,1]        [,2]        [,3]         [,4]
##  [1,]  1.01323287 -0.91056130  0.22300841  0.978403744
##  [2,]  0.27599802 -0.37879572 -1.01693690  0.762334090
##  [3,] -0.03203048 -0.51511158  0.79037437  0.002639731
##  [4,] -0.09193560  0.07907323  0.14790265 -0.039062606
##  [5,] -0.73131494  0.07549138  0.17234648 -0.009426905
##  [6,] -0.50890530  0.18561143  0.72008751  0.496708973
##  [7,] -0.01427086 -0.28208810 -0.11196553 -0.537226514
##  [8,]  0.25274376  0.15911871 -0.11647583  0.011376558
##  [9,] -0.21570241 -0.33741651  0.47711231 -0.459624091
## [10,] -0.23138642 -0.92727329  0.04955606 -0.017246277
## [11,] -0.19430907 -0.38257361 -0.26531432  0.630049848
## [12,] -0.66903676  0.17674572 -0.16014889  0.067177474
## 
## [[4]]
## , , 1
## 
##            [,1]        [,2]        [,3]      [,4]
## [1,]  0.3666743  0.86477950  0.61013519 0.4585473
## [2,] -0.7287975  0.89967854 -0.26984724 0.1715412
## [3,] -0.2753835 -0.27174429 -0.02783652 0.5228947
## [4,]  0.2106197 -0.09700542 -0.95923363 0.2126882
## 
## , , 2
## 
##            [,1]        [,2]       [,3]        [,4]
## [1,] -0.5464374 -0.57925350  0.2394385  0.13922962
## [2,] -0.2060458  0.01186567 -0.2479210  0.52160749
## [3,]  0.6099824 -0.40519047 -0.1605294 -0.05806732
## [4,]  0.7619924 -0.03927812  1.0184070  0.76831379
## 
## , , 3
## 
##            [,1]       [,2]       [,3]          [,4]
## [1,] 0.01579773 -0.0176360  0.1042295 -0.6998017767
## [2,] 1.12913619 -0.1744167  0.1694019 -0.0537070606
## [3,] 0.11367376 -0.5660887  0.5790242 -0.0007441308
## [4,] 0.00766533 -0.7063434 -0.6071945 -0.0742630499
## 
## , , 4
## 
##            [,1]        [,2]        [,3]       [,4]
## [1,] -0.4415180 -0.45882179 0.189424993 -0.4583560
## [2,]  0.4881191 -0.07763796 0.534154371 -0.3413463
## [3,] -0.1194739  0.71146855 0.622996667  0.7601848
## [4,] -0.2681222  0.75851682 0.004319161 -0.7493150
## 
## , , 5
## 
##            [,1]        [,2]       [,3]      [,4]
## [1,]  0.3594965 -0.44704612  0.1057529 0.6412851
## [2,] -0.4734836 -0.02154132  0.8362356 0.0703729
## [3,] -0.1015575  0.59180270 -0.2234910 0.3626638
## [4,]  0.7883569  0.48251087  0.6203973 0.1364310
## 
## , , 6
## 
##            [,1]       [,2]       [,3]       [,4]
## [1,]  0.6413329 -0.0520821  0.0313601 -0.7862978
## [2,]  0.2766572  0.5259953  1.0408831  0.7330976
## [3,] -0.4612240  0.6870955 -0.2413829 -0.1735197
## [4,] -0.8082452 -0.7595234  0.7075496 -0.6229779
## 
## 
## [[5]]
##              [,1]        [,2]        [,3]        [,4]
##  [1,] -0.11900916  0.72971816 -0.33534413 -0.16856932
##  [2,] -0.56655781 -0.01421153 -0.44314434  0.25331983
##  [3,]  0.02974235  0.39626672  0.35083549 -0.16975625
##  [4,]  0.08175985  0.01612911  0.06762713  0.03284772
##  [5,] -0.10161875 -0.04482800  0.34340020  0.33128165
##  [6,]  0.27189293  0.34943436  0.27667359  0.51511528
##  [7,] -0.22418646 -0.17966680  0.14259016 -0.36500830
##  [8,]  0.11894101 -0.03845361 -0.18004009 -0.05569821
##  [9,] -0.06241772  0.05742287  0.38727606 -0.27966268
## [10,] -0.59027383  0.15893478  0.33824411 -0.16558486
## [11,] -0.38817067  0.16736336 -0.01299629  0.34651855
## [12,] -0.16592813 -0.15620554  0.18829157  0.39244612
## 
## [[6]]
##             [,1]         [,2]         [,3]         [,4]
## [1,]  1.06913329 -0.072243449  0.064732772 -0.056487174
## [2,] -0.07224345  0.874066557  0.007063041 -0.007362562
## [3,]  0.06473277  0.007063041  0.948800593 -0.009399430
## [4,] -0.05648717 -0.007362562 -0.009399430  1.001389459
## 
## [[7]]
##   [1]  1.0132328710 -0.9105612976  0.2230084069  0.9784037440  0.2759980196
##   [6] -0.3787957167 -1.0169369015  0.7623340902 -0.0320304800 -0.5151115759
##  [11]  0.7903743678  0.0026397307 -0.0919356017  0.0790732281  0.1479026489
##  [16] -0.0390626056 -0.7313149438  0.0754913824  0.1723464850 -0.0094269050
##  [21] -0.5089053016  0.1856114264  0.7200875059  0.4967089726 -0.0142708566
##  [26] -0.2820881047 -0.1119655280 -0.5372265136  0.2527437636  0.1591187075
##  [31] -0.1164758339  0.0113765575 -0.2157024076 -0.3374165063  0.4771123077
##  [36] -0.4596240914 -0.2313864218 -0.9272732890  0.0495560606 -0.0172462772
##  [41] -0.1943090688 -0.3825736088 -0.2653143217  0.6300498485 -0.6690367591
##  [46]  0.1767457169 -0.1601488935  0.0671774739  0.3666742930  0.8647795015
##  [51]  0.6101351924  0.4585473120 -0.7287975411  0.8996785444 -0.2698472400
##  [56]  0.1715411534 -0.2753834548 -0.2717442880 -0.0278365186  0.5228947384
##  [61]  0.2106197266 -0.0970054220 -0.9592336331  0.2126882478 -0.5464374357
##  [66] -0.5792534966  0.2394384985  0.1392296218 -0.2060458352  0.0118656721
##  [71] -0.2479210148  0.5216074862  0.6099824149 -0.4051904672 -0.1605293560
##  [76] -0.0580673216  0.7619923905 -0.0392781166  1.0184070112  0.7683137886
##  [81]  0.0157977321 -0.0176359999  0.1042294944 -0.6998017767  1.1291361885
##  [86] -0.1744166518  0.1694018671 -0.0537070606  0.1136737596 -0.5660887218
##  [91]  0.5790242132 -0.0007441308  0.0076653298 -0.7063434037 -0.6071945083
##  [96] -0.0742630499 -0.4415179511 -0.4588217919  0.1894249926 -0.4583559714
## [101]  0.4881190547 -0.0776379596  0.5341543712 -0.3413463261 -0.1194739358
## [106]  0.7114685497  0.6229966674  0.7601847997 -0.2681222285  0.7585168195
## [111]  0.0043191611 -0.7493149802  0.3594964704 -0.4470461217  0.1057529089
## [116]  0.6412850683 -0.4734836453 -0.0215413243  0.8362356329  0.0703728975
## [121] -0.1015575348  0.5918026977 -0.2234910091  0.3626638219  0.7883569250
## [126]  0.4825108655  0.6203972867  0.1364309830  0.6413329012 -0.0520821005
## [131]  0.0313601037 -0.7862978038  0.2766571825  0.5259952870  1.0408830938
## [136]  0.7330976036 -0.4612239907  0.6870955302 -0.2413828743 -0.1735197274
## [141] -0.8082452080 -0.7595233690  0.7075495824 -0.6229778721
## 
## [[8]]
## [1] -5.625305

Note that the parameter values for mining the matrix were identical to those chosen in the simulations from 19 AUG 2021. The results are as expected; all traits evolve in response to selection to feed on both crop types and resist both pesticide types.

Asexual, no rotation