Fixed empty runway bug

I’ve fixed the error thrown when trying to evalueate the furthest occupied pot in an extinct replicate. This uses a new function, last_occupied_pot(), in helpers.R, which is called from within the model function but can also be called from the analysis code (see example usage in make_Ler_sims.R). The optional argument zero allows the return value for an empty replicate to be something other than zero (to avoid the risk of zero-index values).

I ran make_Ler_sims successfully, and the warnings indicated that it ran into 26 instances of an extinct replicate.

Code profiling

I’ve profiled the model using make_Ler_sims. The results are in logs/make_Ler_sims-2022-10-03.Rprofvis (not on git as it’s large!)

• Total time in iterage_genotype(): 139300 ms
• The bulk of that is in seed_sampling(): 134360
• Within that, the bulk is in the aaply call to rgengamma (25k), and the aaplys to disp_table (50k each)
• Within dist_table(), 75k is in the line raw_table <- table(dists). This is more than half the total time!

So the first thing to do is to find a faster replacement for table. This StackOverflow page suggests that tabulate() (in base R) is about 14 times faster than table(), which would bring this operation down to 5-6k ms. According to the help page, “tabulate is the workhorse for the table function.” So the time is spent in all the robustness of table. It may take some work to figure out exactly how to call tabulate.

Note that this is not one of the issues I flagged in my manual code review!

Code profiling

I’ve changed disp_table() to use tabulate rather than table, and confirmed that they return the same results. Microbenchmarking (see test_disp_table()) suggests that there is a 100-fold speedup in that function!

Profiling make_Ler_sims.R with this change now speeds things up a lot (make_Ler_sims-2022-10-04.Rprofvis). About half of the total time is in calculating the gengamma random numbers, which doesn’t make use of the internal vectorizing, instead using aaply to go across parameters. Table building now takes about 40% of total time, but less than 10% of that seems to be actually inside disp_table; apparently aaply has a huge performance overhead! It seems that the base apply may be much faster, but I wonder whether the bugs I was trapping by using .drop = FALSE will reappear. It might also be that those commented-out aperm statements were applied to an earlier attempt at using apply.

General

I updated my home mac to R 4.4.1, and updated all libraries.

I found the code to revert blogdown posts to their old structure, and added it to the project’s .Rprofile. I also added the option not to automatically rebuild the website every time a file is saved, as that is problematic if there are code errors and also takes nontrivial time with the (uncached) ProjectTemplate loading.

Bugs

I removed the empirical = TRUE arguments to mvrnorm().

Efficiency

Let’s look at replacing aaply with apply. Things to check include the shape of the returned function and what happens to dimensions with value 1. The output uses up lots of space so I don’t show it; the comments below summarize the results of running the code.

xx <- array(1:24, dim = c(3,4,2))
# Scalar function return
y1 <- aaply(xx, c(1,2), sum)
y2 <- apply(xx, c(1,2), sum)
y1
y2
class(y1)
class(y2)
# Vector function return
y1 <- aaply(xx, c(1,2), function(x) c(sum(x), prod(x)))
y2 <- apply(xx, c(1,2), function(x) c(sum(x), prod(x)))
y1
y2
aperm(y2, c(2, 3, 1))
class(y1)
class(y2)

xx <- array(1:8, dim = c(1,4,2))
xx
aaply(xx, c(1,2), function(x) c(sum(x), prod(x)), .drop = FALSE)
apply(xx, c(1,2), function(x) c(sum(x), prod(x))) %>% aperm(c(2, 3, 1))
aaply(xx, c(1,2), function(x) c(sum(x), prod(x)))

xx <- array(1:8, dim = c(4,1,2))
xx
aaply(xx, c(1,2), function(x) c(sum(x), prod(x)), .drop = FALSE)
apply(xx, c(1,2), function(x) c(sum(x), prod(x))) %>% aperm(c(2, 3, 1))
aaply(xx, c(1,2), function(x) c(sum(x), prod(x)))

xx <- array(1:8, dim = c(2,4,1))
xx
aaply(xx, c(1,2), function(x) c(sum(x), prod(x)), .drop = FALSE)
apply(xx, c(1,2), function(x) c(sum(x), prod(x))) %>% aperm(c(2, 3, 1))
aaply(xx, c(1,2), function(x) c(sum(x), prod(x)))

• When the function returns a scalar, the outcomes are the same, except that rows and columns from apply are unnamed. The latter should not be an issue for me.
• When the function returns a vector, that vector is put into the leading dimension of the returned array by apply, whereas it remains in its original position under aaply. This is clearly why I had those aperm statements.
• When one of the dimensions being parsed has value 1, then aaply with the defaults drops a dimension. But apply by itself seems fine.

So it seems that apply %>% aperm should be a drop-in replacement for aaply(.drop = FALSE). So why did I make the switch in the first place? Was I just trying to be modern, or does aaply have some other advantage? I made the change on 5/29/19, in a commit that has massive numbers of code updates, with no comments on individual changes. It may just be that I liked the elegance of skipping the aperm step.

I’ve tried adding the apply version back in (I had to fix the order of indices in aperm); compared with aaply, it’s about 2.5 times faster. It’s still the case that most of that time is in the overhead. apply has a lot of cases, but if I’ve parsed it correctly, in my use case it’s doing a for loop over the nreps x npots cases. Further, rather than calling the function directly in that loop, it’s calling something called forceAndCall which is trapping for closures, which may add more time. And it’s hard to know how much is being added by aperm. This post suggests that apply is as fast as a for loop, but the system times were very crude (0.02 sec vs. 0.00 sec). I run a microbenchmark based on their example:

obs <- array(rnorm(10*13*12), c(10,13,12))
floop <- function(obs) {
  l2 = dim(obs)[2] 
  l3 = dim(obs)[3]
  minObs<-matrix(nrow=dim(obs)[2],ncol=dim(obs)[3])
  for (i in 1:l2) 
    for (j in 1:l3) 
      minObs[i,j]<-min(obs[,i,j],na.rm = TRUE)
  return(minObs)
}
mbm <- microbenchmark(aaply(obs, c(2,3), min), floop(obs), apply(obs, c(2,3), min))

Unit: microseconds
                     expr       min        lq       mean     median        uq       max neval
 aaply(obs, c(2, 3), min) 26991.818 27617.330 29557.6689 29839.2335 30930.627 35212.618   100
               floop(obs)   219.243   245.001   285.7647   264.8650   303.561   439.593   100
 apply(obs, c(2, 3), min)   281.425   317.083   353.9656   341.7015   372.688   573.770   100

In this case, there’s only a 20% speedup of the for loop vs apply. So maybe it’s in the aperm call?

mbm2 <- microbenchmark(apply(obs, c(2,3), min), apply(obs, c(2,3), min) %>% aperm(c(2,1)))

Unit: microseconds
                                        expr     min       lq     mean  median       uq     max neval
                    apply(obs, c(2, 3), min) 273.826 288.7125 324.2954 299.926 339.1565 482.853   100
 apply(obs, c(2, 3), min) %>% aperm(c(2, 1)) 293.200 310.4655 345.1462 324.005 343.2050 554.648   100

Here there’s basically no difference. So the overhead is all in the for loop, and we probably can’t get away from it without working out a direct vectorization. But it is curious that the speedup from aaply to apply in the microbenchmark is so much greater than it is in the profile.

At any rate, I finalize the code with apply and aperm.

Efficiency

I just profiled apply vs aaply for the random number generation, and the former takes about half as long. So I’ll go with that.

Next, I compared a single call to rgengamma with a long vector of parameters to the apply call. Even with for loops before and after to simulate preparing and parsing the data, there was more than a 10-fold speedup (see test_rgengamma()). So it seems worthwhile to pursue this.

Haven’t quite gotten there, but here are some code snippets that point where I’m heading:

disp_seeds <- matrix(c(5, 4, 0, 2, 10, 3), nrow = 2, byrow = TRUE)
QQ <- matrix(1:6, nrow = 2, byrow = TRUE)
x <- array(c(disp_seeds, QQ), dim = c(2,3,2))
y <- unlist(apply(x, c(2,1), function(x) rep(x[2], x[1]), simplify = FALSE))
max_ds <- max(disp_seeds) # to set the array dimension to pad to
yy <- array(0, dim = c(2,3,max_ds))

Then a double for loop to fill in appropriate parts of yy with appropriate parts of y. Note that the apply has the dimensions reversed to ensure that y comes out row-wise.

Efficiency

OK, I’ve implemented the revised code. The savings aren’t as dramatic as microbenchmark had suggested (perhaps because of memory issues?); maybe a 2-fold speedup. But it’s something!

Code inconsistencies

I’ve rationalized the way dispersal tails are handled by making det_kernel simply calculate out to the value of controls$max_pots and getting rid of the reflection in seed_sampling. We still need a rationalized value of maximum dispersal distance.

Conceptual issues

I’ve added a value to pass into iterate_genotype to use for the temporal component of environmental stochasticity in seed production, with a sensible default. This then gets passed through to ES_seeds. This runs with both test_iterate_genotype and with an appropriately updated make_Ler_sims

Documentation

• Cleaned up comments
• Fixed equation derivation in gompertz_seeds()
• Added documentation for missing input list values

Code structure

I have separated the model inputs into plant_params, expt_params, and sim_settings. Everything is appropriately propagated through model.R, and I have updated test_iterate_genotype() to work with the new structure. I’ll need to update some of the other test functions, as well as the simulation drivers.

I did realize on doing this that I’m often passing the entire list through simply to get a single value; in those cases I should simplify things.

• I’ve pulled n_pots out of the control arrays, as it is only used internally and is updated dynamically.
• I’ve simplified the call to gapify()
• Simplified call to seed_sampling()
• Simplified call to det_kernel()
• Simplified call to kernel_stoch()
• Simplified call to DS_seeds()
• Simplified call to ES_seeds(), and fixed mistake in application of ES_seed_time
• Simplified call to Gompertz_seeds()
• That’s everything for eliminating unessential list-passing!

New bugs

• There seems to be an issue with apply() dropping an array dimension if there is only one replicate that hasn’t gone extinct.

New bugs

In trying to run the visualization code, I was getting errors about on out-of-date grouping scheme for LerC_spread. I tried re-munging, but that threw an error in 04-Ler_dispersal_stats.R. So apparently the package updated have left something non-functional in the munge code that I’ll need to update!