Rewrite the exponential decay process in C++.

[plietar]

The existing process written in R needs to copy the contents of each variable from C++ using v$get_values(), then after scaling the vector it would copy the result back into C++ using v$queue_update.

The amount of data copied and the time it took was pretty significant. 6 double variables, one for each kind of immunity, need to be updated in full at each time step, each as big as the population size.

Moving this into C++ removes the need for any copy at all, besides the multication loop. Values are read out of a reference to the vector held by the DoubleVariable, the result of the multication is moved to the queue, and finally individual moves the vector in the queue into the DoubleVariable. The speedup from this change for a 1M population size is around 10%.

An alternative optimization I considered was to compute the exponential decay lazily, recording only the timestep and value at which the immunity was last updated and using the closed form expression of the exponential decay. This would avoid the need to have mass updates of the immunity variables at every time step. Unfortunately in my testing this ends up being slower than even the current implementation, with all the time being spent in calculating the current value. This would also be a much more intrusive change, since every use of the immunity variables needs to be modified to take the last update timestep, the current timestep and the decay rate into consideration.

[plietar]

The broken tests are because of the new individual version, will be fixed by mrc-ide/individual#192

[giovannic]

This is great stuff!

Though not necessary, it would be good to see the performance improvement compared to a naive copy, modify, move implementation (which wouldn't require all the iterator code).

Another alternative to potentially consider: Rounding off low immunities to zero and excluding them from the update? Milage would likely depend on the threshold for low immunity and the transmission level.

[giovannic]

What about copy, modify, move?

[giovannic]

It's very surprising that gcc doesn't optimise the second approach. Is that because of the default optimisation level set by R packaging? Or is this generally the best way to do it?

[plietar]

Is that because of the default optimisation level set by R packaging?

Not really, I get similar results even at -O3 -march=alderlake.

See the disassembly at https://godbolt.org/z/9jWKbTev7. There's a lot going on, but even just at a glance the copy_modify and zero_initialized functions each start with a call to memmove or memset, respectively. reserve_pushback doesn't have any such call, but it doesn't have any unrolling and vectorization. The only multiplication is a single vmulsd instruction, which only multiplies one double value per call, whereas the other functions use multiple calls to vmulpd ymm, which does four multiplications at once (along with some vmulpd xmm and vmulsd calls to round things off).

Clang is even more aggressive about unrolling, and generates 4 vmulpd ymm instructions for loop iteration, and generally produces an easier to read assembly. Still doesn't optimize out the memmove and memset calls though.

Here are the microbenchmarks cited in the comment https://gist.github.com/plietar/72fc5160f79e19379418ce391db1624a

[plietar]

To be honest I am a little on the fence about whether to use the iterator stuff or not. It was about 2% speedup on the overall simulation compared to the zero-initialization code, which is small but not insignificant. Given that all the complexity is contained within that one function and doesn't leak out, I think I'm inclined to keep it.