piLaboratory / GillesCom

Gillespie Simulation of Ecological Communities Dynamics
2 stars 0 forks source link

Transfer bdm code to C++ #1

Closed andrechalom closed 8 years ago

andrechalom commented 8 years ago

Transferred bdm code to C++. Temporarily renamed bdm function to Rbdm in order to run benchmarks. C implementation is 10 times faster:

> Rf <- function(N0) {for (i in 1:1000) N0 <- Rbdm(N0, alpha)} # R implementation
> Cf <- function () {for(i in 1:1000) bdm()} # C++ implementation
>
> compare <- function(N) {
+ Interaction(N)->alpha
+ N0 <- rep(1, N)
+ Init_Community(N0, alpha)
+ microbenchmark(Rf(N0), Cf())
+ }
> compare(100)
Unit: milliseconds
 expr       min        lq      mean   median        uq       max neval cld
 Rf() 106.85142 114.49811 122.59260 117.6926 120.55049 228.19975   100   b
 Cf()  14.31535  14.43105  14.90403  14.5538  14.67543  16.59435   100  a 
> compare(10)
Warning: overwriting previous Community
Unit: milliseconds
   expr       min        lq      mean    median        uq       max neval cld
 Rf(N0) 21.481973 22.209533 23.398928 22.834599 24.233413 28.677935   100   b
   Cf()  2.632542  2.740774  3.153657  2.808976  3.001187  6.668018   100  a 

Now we need to run some optimizations on the C code.

andrechalom commented 8 years ago

New improvement, now bdm accepts a "count" argument, so the loop is made inside C code:

> microbenchmark(for(i in 1:50000) bdm(), bdm(50000))
Unit: milliseconds
                     expr      min       lq     mean   median       uq      max neval cld
 for (i in 1:50000) bdm() 678.1162 688.9477 692.2319 690.8555 692.5461 793.9037   100   b
               bdm(50000) 556.0742 557.7395 558.3793 558.2241 558.7047 561.8939   100  a 
andrechalom commented 8 years ago

"Caching" the death rate slope (b-d0)/K on the constructor instead of repeating this every time gave me more 5% speedup. Now I don't see how can we improve it more, as the N vector needs to be recalculated at every step (because changing any abundance changes every position in N), and the multiplication step in abundance * interaction is responsible for over 90% of the running time.

Update: more 2% speedup as I collapsed two lines of code, now d is calculated directly, without needing N.

Waiting checks from PI to close this and remove Rbdm() from source.

andrechalom commented 8 years ago

Removed from main branch; latest version with original bdm and run.bdm is 24a50ad