Benchmarks of various scientific programming languages for theoretical modelling

For a theoretical modelling course for biology students, I am trying to decide which would be the best technical programming language for doing evolutionary simulations in terms of elegance and compactness, whilst still being performant (this is for simulations, for analytical stuff I am using Mathematica).

As a simple example I coded a simulation module that looks at the evolution of a continuous trait in an asexually reproducing population under density dependent competition in discrete time (i.e. using non-overlapping generations, using recurrence equations, where parents all reproduce simultaneously, produce offspring, after which the parents all die and the offspring become new parents) in 4 different languages: R, Julia, Mathematica and Matlab.

None, however, appear very fast, and none scale up to, say, simulating 1 million generations in my simple implementation. Any thoughts on optimising code performance, memory usage or scalability in either implementation would be welcome, as well as any recommendations or ports to other high-level technical computing languages like Python/Numba/Cython. An elegant C, C++, Java or JVM OO-based programming language implementation would be cool, too.

R implementation:

```
# function to calculate new offspring trait values from parent population trait values tv
# given a density-dependent fitness function fitnessfunc,
# a mutation rate mutrate and standard deviation of mutational effects stdev
doStep = function (tv, fitnessfunc, mutrate, stdev) {
n = length(tv) # current population size
numberoffspring = rpois(n, fitnessfunc(R=tv, popsize=n)) # number of offspring of each parent
newtv = rep(tv, times=numberoffspring) # offspring are copies of their parents
# save occasional mutation which we apply below
n = length(newtv)
nrmutants = rbinom(1, n, mutrate)
rnoise = rnorm(nrmutants, mean=0, sd=stdev)
rndelem = sample(1:n, nrmutants, replace=FALSE)
newtv[rndelem] = newtv[rndelem] + rnoise
return(pmax(newtv

I started with the Julia code you had and also got ~20 seconds, so I think my timings are similar to yours. Let me give a step by step breakdown on how to do this.

To start, notice that if you are running code in the REPL that variables defined there are global. This incurs a good performance cost. There are two ways to deal with this:

1) Wrap it all in a function like:

function plotThePops()
   #all code in here
end
plotThePops()

or 2) Change some globals to constants.

The compiler needs this because otherwise it cannot properly infer the types of these variables (at global scope, they can change anywhere!) which means that it needs to compile a bunch of helper code everywhere in case the types change.

I went with the const route, to get the following code:

using Distributions
# Pkg.clone("https://github.com/ChrisRackauckas/VectorizedRoutines.jl") # for R.rep function
# Will be Pkg.add("VectorizedRoutines") after being added to package system
using VectorizedRoutines

function doStep(tv, fitnessfunc, mutrate, stdev)
    n = length(tv) # current population size
    numberoffspring = [rand(Poisson(theta)) for theta in fitnessfunc(tv, n)] # number of offspring of each parent
    newtv = R.rep(tv, each = numberoffspring) # offspring are copies of their parents
    # but some of them mutate, so we also apply mutation
    n = length(newtv)
    nrmutants = rand(Binomial(n, mutrate),1)[1]
    rnoise = randn(nrmutants)*stdev
    rndelem = sample(1:n, nrmutants[1], replace=false)
    newtv[rndelem] = newtv[rndelem] + rnoise
  return(max(newtv,0))
end

function evolve(init_traitv, fitnessfunc, mutrate, stdev, ngens)
    Results = Vector{Vector{Float64}}(ngens)
    Results[1] = init_traitv
    for idx = 2:ngens
        Results[idx] = doStep(Results[idx - 1], fitnessfunc, mutrate, stdev)
    end
    return(Results)
end

const psize = 1000 # population size;
const ngens = 10000  # nr of generations to simulate;
const mutrate = 0.005 # mutation rate;
const stdev = 0.05 # st dev of mutant trait values;
const k = 2*psize # carrying capacity;
srand(1);
const init_traitv = rand(2.5:.1/psize:2.6,psize) # initial trait values;
fitnessfunc(R, popsize; K=k) = max((1 + R*(1 - popsize/K)), 0)

@profile res = evolve(init_traitv, fitnessfunc, mutrate, stdev, ngens)
@time res = evolve(init_traitv, fitnessfunc, mutrate, stdev, ngens)

popsizes = Int[length(res[idx]) for idx = 1:length(res)] # population size;
densities = hcat([hist(res[idx],0:3/100:3)[2]/length(res[idx]) for idx = 1:length(res)]...)' # densities;

## Pkg.add("Plots")
## Pkg.add("GR") if you want to use the GR plotting library.
using Plots
gr()
p1 = plot(popsizes, 1:ngens, lab = "",
          xlabel = "Population size", yaxis = ("Generation",(0,ngens)));
p2 = heatmap(densities, xticks = (100/6:100/6:100, 0.5:0.5:3), fill = :grays,  # axis labels could probably be specified more elegantly
             xlabel = "Reproductive rate R");
plot(p1, p2)
gui()

Note that I also switched to the GR plotting backend which focuses on speed. That's not in the timings, but just thought you should know. This reduces the timings down to around ~17 seconds. Not much, but okay.

So where do we go next? To find out, change the timing line to profiling:

@profile res = evolve(init_traitv, fitnessfunc, mutrate, stdev, ngens)
using ProfileView
ProfileView.view()

This opens up a window like this:

Here's how you read that diagram. The block length denotes the percentage of the time which is taken up by a line, and if you scroll over the block it tells you which line it is. If a block is on top of another, that means the code on top was called from the code below it. The line that I circled is for the dostep function, and the part that I circled was for the following line:

numberoffspring = [rand(Poisson(theta)) for theta in fitnessfunc(tv, n)]

This means that ~80% of the time in your timing was due to this line! Obviously that means it was implemented inefficiently. Indeed, this can be improved by unrolling and turning off bounds checking. So I implemented the following function in VectorizedRoutines.jl:

function rpois(n::Int,p::Vector{Float64})
  out = Vector{Int}(n)
  for i in 1:n
    @inbounds out[i] = rand(Poisson(p[i]))
  end
  out
end

and then the offending dostep line was changed to:

numberoffspring = R.rpois(n,fitnessfunc(tv, n))

(if you want put the function in your script you could just call it as rpois). Now the times are down to 2.8-4 seconds depending on randomness (and ~32-34 seconds for ngens = 100000). Profiling again I get

That big red block is for VectorizedRoutines.jl's R.rep function. Where it's caught up is:

for j in eachindex(x)
    @inbounds for k in 1:each[j]
      index += 1
      @inbounds rval[index] = x[j]
    end
  end

This is just straight moving numbers around and thus probably can't be improved.

Thus the final code is this:

```
using Distributions
# Pkg.clone("https://github.com/ChrisR

StackExchange Code Review Q#134926, answer score: 21