# 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,0))
}

# function to iterate this a given number of generations, starting with trait values tv
evolve = function (tv, fitnessfunc, mutrate, stdev, ngens) {
sapply(1:ngens, function (gen) { tv <<- doStep(tv, fitnessfunc, mutrate, stdev) # slightly dodgy solution by defining tv globally, but works
return(tv) })
}

# example parameters
psize = 1000 # pop size
ngens = 10000  # nr of generations to simulate
mutrate = 0.005 # mutation rate
stdev = 0.05 # st dev of mutant trait values
k = 2*psize # carrying capacity
init_traitv = runif(psize, 2.5, 2.6) # initital trait values
fitnessf = function (R, popsize, K=k) pmax((1 + R*(1 - popsize/K)), 0) # density-dependent fitness function

# do simulation and plot results
set.seed(1)
system.time(out <- evolve(init_traitv, fitnessf, mutrate, stdev, ngens))

popsizes = sapply(out,function(x) length(x))
# densities = sapply(out, function (x) hist(x, breaks=seq(0,3,length.out=100), plot=FALSE)\$count/sum(x,na.rm=TRUE) )
densities = sapply(out, function (x) tabulate(findInterval(x, vec=seq(0,3,length.out=100)),nbins=99)/sum(x,na.rm=TRUE) ) # slightly faster
par(mfrow=c(1,2))
plot(popsizes, 1:ngens, type="l", xlab="Population size", ylab="Generation", col="steelblue", las=1, cex.axis=0.7, ylim=c(1,ngens), yaxs="i")
image(x=seq(0,3,length.out=99), y=1:ngens, z=densities, col=colorRampPalette(c("grey100","grey0"))(50),
useRaster=TRUE, xlim=c(0,3), las=1, cex.axis=0.7, xlab="Reproductive rates (R)", ylab="Generation")
box() Julia implementation:

## Pkg.add("Distribution")
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, k)] # 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)
rnoise = randn(nrmutants)*stdev
rndelem = sample(1:n, nrmutants, replace=false)
newtv[rndelem] = newtv[rndelem] + rnoise
return(max(newtv,0))
end

function evolve(tv, fitnessfunc, mutrate, stdev, ngens)
Results = Array(Array, ngens)
Results = tv
for idx = 2:ngens
Results[idx] = doStep(Results[idx - 1], fitnessfunc, mutrate, stdev)
end
return(Results)
end

psize = 1000 # population size;
ngens = 10000  # nr of generations to simulate;
mutrate = 0.005 # mutation rate;
stdev = 0.05 # st dev of mutant trait values;
k = 2*psize # carrying capacity;
srand(1);
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);

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

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

using Plots
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) Mathematica implementation:

(* fast Poisson random number generator *)
fastPoisson = Compile[{{lambda, _Real}},
Module[{b = 1., i, a = Exp[-lambda]},
For[i = 0, b >= a, i++, b *= RandomReal[]]; i - 1],
RuntimeAttributes -> {Listable}, Parallelization -> True];

(* function to mutate trait values *)
mutate[traitvalues_, mutrate_, stdev_] :=
Module[{tv, n, nrmutants, rnoise, rndelem},
tv = traitvalues;
n = Length[tv];
nrmutants = RandomVariate[BinomialDistribution[n, mutrate]];
rnoise = RandomReal[NormalDistribution[0, stdev], nrmutants];
rndelem = RandomChoice[Range[n], nrmutants];
tv[[rndelem]] += rnoise;
Clip[tv, {0, 10^100}]]

(* function to calculate new offspring trait values from parent \
population with trait values *)
doStep[tv_, fitnessfunc_, mutrate_, stdev_] :=
Module[{n, fitnessinds, numberoffspring, newtv},
n = Length@tv;
fitnessinds = fitnessfunc[tv, n];
numberoffspring = fastPoisson[fitnessinds];
newtv = Join @@ MapThread[ConstantArray, {tv, numberoffspring}];
mutate[newtv, mutrate, stdev]]

(* function to iterate a number of generations *)
evolve[fitnessfunc_, initpop_, mutrate_, stdev_, generations_] :=
NestList[doStep[#, fitnessfunc, mutrate, stdev] &, initpop,
generations]

(* function to plot results of one run *)
PlotResult[traitvalues_, maxphen_] :=
Module[{generations, pop, maxscaleN, minscaleN, maxscaleP,
frequencydata},
generations = Length[traitvalues];
pop = Length /@ traitvalues;
maxscaleN = Max[pop];
minscaleN = Min[pop];
maxscaleP = Max[Max[traitvalues], maxphen];
frequencydata = (BinCounts[#, {0, maxscaleP, 1/100}] & /@
traitvalues)/(pop + 0.00001);

GraphicsRow[{ ListPlot[Transpose[{pop, Range[generations]}],
Joined -> True, Frame -> True,
FrameLabel -> {"Population size", "Generation"}, AspectRatio -> 2,
PlotRange -> {{Clip[minscaleN - 50, {0, \[Infinity]}],
maxscaleN + 50}, {0, generations}}],
Show[ArrayPlot[frequencydata,
DataRange -> {{0, maxscaleP}, {1, generations}},
DataReversed -> True], Frame -> True, FrameTicks -> True,
AspectRatio -> 2,
FrameLabel -> {"Phenotype frequency", "Generation"}] }]
]

(* example parameters and simulation *)
psize = 1000; ngens = 10000; mutrate = 0.005; stdev = 0.05; k =
2*psize;
f[R_, popsize_] =
Clip[(1 + R (1 - popsize/k)), {0.000001, 10^100}];
First@AbsoluteTiming[
traitvalues =
EvolveHapl[f, RandomReal[{2.5, 2.6}, psize], mutrate, stdev,
ngens];]
PlotResult[traitvalues, 3] Matlab implementation:

doStep.m function file:

function [ tv ] = doStep( tv, fitnessf, mutrate, stdev )
% Function to calculate new offspring trait values from parent population trait values tv
% given a density-dependent fitness function fitnessf,
% a mutation rate mutrate and standard deviation of mutational effects stdev

n = length(tv); % current population size
numberoffspring = poissrnd(fitnessf(tv, n)); % number of offspring of each parent
newtv = repelem(tv, numberoffspring); % offspring are copies of their parents
% save occasional mutation which we apply below
n = length(newtv);
nrmutants = binornd(n, mutrate);
rnoise = normrnd(0, stdev, 1, nrmutants);
rndelem = datasample([1:n], nrmutants, 'Replace', false);
newtv([rndelem]) = newtv([rndelem]) + rnoise;
tv = max(newtv,0);

end


evolve.m function file:

function [ Results ] = evolve( tv, fitnessf, mutrate, stdev, ngens )
% Iterate evolution given number of generations and output all trait
% vectors in each generation as cell array Results
out{1} = tv;
for idx = 2:ngens
out{idx} = doStep(out{idx - 1}, fitnessf, mutrate, stdev);
end
Results = out;

end


Matlab script file:

% example parameters
psize = 1000; % pop size
ngens = 10000;  % nr of generations to simulate
mutrate = 0.005; % mutation rate
stdev = 0.05; % st dev of mutant trait values
k = 2*psize; % carrying capacity
rng('default');
init_traitv = 2.5 + (2.6-2.5).*rand(1,psize); % initital trait values
fitnessf = @(R, popsize) max((1 + R.*(1 - popsize/k)), 0) % density-dependent fitness function

% do simulation and plot results
tic
out = evolve(init_traitv, fitnessf, mutrate, stdev, ngens);
toc

popsizes = cellfun(@(el) length(el), out);
densities = reshape(cell2mat(cellfun(@(el) histcounts(el, 0:3/100:3), out, 'UniformOutput', false)),100,ngens)';

% plot results
clf
subplot(1,2,1)
plot(popsizes,1:ngens)
xlabel('Population size', 'FontSize', 16)
ylabel('Generation', 'FontSize', 16)
subplot(1,2,2)
imagesc(linspace(0,3,100), linspace(1,ngens,100), densities)
set(gca,'ydir','normal')
colormap(flipud(gray))
xlabel('Reproductive rate (R)', 'FontSize', 16)
% colorbar()
% ylabel(p1,'Density','FontSize', 16) Timings for these 4 implementations on my Intel i7 2.4 GHz 8 Gb laptop under Windows 8.1:

With ngens = 10000 and psize = 1000:

• R 3.2.3: 5.8 s (Microsoft R Open, with multithreaded BLAS/LAPACK libraries and 4 cores but not sure they're used)
• Julia 0.4.6: 20.6 s
• Matlab 2015a: 31.9 s
• Mathematica 10: 96.5 s

With ngens = 100000 and psize = 1000:

• R 3.2.3: 38 s (Microsoft R Open, with multithreaded BLAS/LAPACK libraries and 4 cores but not sure they're used)
• Julia 0.4.6: 231 s
• Matlab 2015a: 328 s
• Mathematica 10: 966 s

From this, R would appear to be the winner, although the relatively poor Julia timings surprised me, given that published benchmarks are generally really good (close to C++), and generally much faster than R.

Any thoughts what I'm doing wrong in my Julia implementation (problems with type stability or due to lack of type declarations)? Note that the plotting routine ideally should be capable of plotting a 1 million generation output, although a bit of downsampling would be acceptable for the heatmap.

• Consider profiling to find out where the potential hotspots are--this may help in focusing on the problem: adv-r.had.co.nz/Profiling.html See also adv-r.had.co.nz/Performance.html
– Matt
Jul 15, 2016 at 10:49
• Thx - one bottleneck I just found appears to be the line "newtv = do.call(c, sapply(1:n,function(i) rep(tv[i], numberoffspring[i])))" - any thoughts to do that faster/better? Jul 15, 2016 at 11:14
• Ha just found it - rep() with times argument :-) Jul 15, 2016 at 11:18
• It sounds like your code doesn't work as intended. Can you clarify this?
– Mast
Jul 17, 2016 at 10:11
• What do you think I'm doing wrong? For me the code seems to work... Jul 17, 2016 at 13:57

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)
rnoise = randn(nrmutants)*stdev
rndelem = sample(1:n, nrmutants, 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 = 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)/length(res[idx]) for idx = 1:length(res)]...)' # densities;

## 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/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 = R.rpois(n,fitnessfunc(tv, n))
newtv = R.rep(tv,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))
rnoise = scale!(randn(nrmutants), stdev)
rndelem = sample(1:n, nrmutants, replace=false)
newtv[rndelem] += rnoise
return(max(newtv,0.0))
end

function evolve(init_traitv, fitnessfunc, mutrate, stdev, ngens)
Results = Vector{Vector{Float64}}(ngens)
Results = 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.0 + r*(1.0 - popsize/K)), 0.0) for r in R]

@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)/length(res[idx]) for idx = 1:length(res)]...)' # densities;

## 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()


Some complexity is hidden because I made a fast version of the Poisson random number generation and added it to the VectorizedRoutines.jl package, but even that function is pretty simple. That's what's nice about Julia: those "low-level utility functions" can be written in Julia so it can be quick/easy to make "fast" versions of things directly in Julia, whereas for R/MATLAB/Python these functions are written in C or Fortran.

I hope this serves a good overview on how to optimize Julia code as well. I may repost this to my blog. I wonder what timings/profile you get.

Edit

With help from Fengyang Wang @TotalVerb via Gitter, we improved the R.rep function in VectorizedRoutines.jl. To get this, just run Pkg.update() and change the R.rep line to:

newtv = R.rep(tv,numberoffspring)


and you should have the improvement. On my computer, the code now runs in 1.1 seconds, so more than a 3x speedup. Given what has been posted, it should easily beat R's timing.

For reference, the new profile looks like this: I circled the RNG time. As mentioned in the comments, from the Google Summer of Code page you can see that there is a project to make a repo RNG.jl with faster psudo-rngs. It looks like it already has a few new ones implemented. So if you run this code at the end of the summer, you may be even faster!

Edit 2

There was a slight change to doStep and fitnessfunc in the final code to reduce allocations. We chopped off a bit more allocations and rpois is using the Rmath library for its RNG. At this point, R and Julia are doing basically the same computations.

R is doing well because it's directly calling C libraries, but from these analyses it seems that an optimized C code is likely only 30% better. So the overhead isn't an issue, but the main hold-ups are the actual algorithms (and the RNG). A new timing will probably put Julia and R at about the same timings. It would be better to test Julia vs R on an algorithm that doesn't vectorize so well!

Edit 3

For completeness I made a version with in-place operations. While the previous version matches R in timing, this one will do slightly better because it reduces allocations in a way you can't do with just vectorized function calls. For example, we can change R.rpois to R.rpois! (! being Julia convention for an in-place function) to re-use the newoffspring vector from the previous round by doing:

  function rpois!(n::Int,p::Vector{Float64},out::Vector{Int})
resize!(out,n)
for i in 1:n
@inbounds out[i] = StatsFuns.RFunctions.poisrand(p[i]) #rand(Poisson(p[i])) #
end
end


Here out is newoffspring. Note that the only difference is that we resize before using it, and now have to pass it around everywhere. For this minor adjustment to the code, we reduce allocation. We can do the same for R.rep to R.rep! by

  function rep!{T1,T2}(x::AbstractVector{T1},rval::AbstractVector{T1},each::AbstractVector{T2} = ones(Int,length(x)))
length_old = length(x)
length_old != length(each) && throw(ArgumentError("If the argument 'each' is not a scalar, it must have the same length as 'x'."))
length_out = sum(each) #times =  1
index = length_out
resize!(rval,length_out)
resize!(x,length_out)
for j in length_old:-1:1
@inbounds @simd for k in 1:each[j]
@inbounds rval[index] = x[j]
index -= 1
end
end
x[1:length_out] = rval[1:length_out]
return(length_out)
end


where now we keep a second vector around for holding the intermediate values (it's required and you cannot write directly into x because of the possibility that a given specimen has 0 offspring. If that wasn't possible you could just write directly into x going backwards and get more of a speed up!).

Note that you can cheat and know that, given the problem, rval will not ever get larger than say 5000, and so you can make a special version of this function where you comment out the line resize!(rval,length_out), pre-allocate rval to a size of 5000, and just never have to deal with re-allocating it. That's slightly "cheating" because the function wouldn't generally work with this (it would error anytime rval needs to resize, so I'm not adding that to the library!), and so the answer would require that you define a rep! function that does this, increasing the complexity of the answer. However, this does chop off a slight bit more of time.

Still, this only saves a small bit of time because the majority of the time is still in the RNG. For reference, the script with the in-place function calls is this:

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,numberoffspring,rval)
n = length(tv) # current population size
R.rpois!(n,fitnessfunc(tv, n),numberoffspring)
n = R.rep!(tv,rval,numberoffspring)
# but some of them mutate, so we also apply mutation
nrmutants = rand(Binomial(n, mutrate))
rnoise = scale!(randn(nrmutants), stdev)
rndelem = sample(1:n, nrmutants, replace=false)
tv[rndelem] += rnoise
return(max(tv,0))
end

function evolve(init_traitv, fitnessfunc, mutrate, stdev, ngens)
Results = Vector{Vector{Float64}}(ngens)
Results = init_traitv
numberoffspring = Vector{Int}(5000)
rval = similar(init_traitv)
for idx = 2:ngens
@inbounds Results[idx] = doStep(Results[idx - 1], fitnessfunc, mutrate, stdev,numberoffspring,rval)
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)

• Comments are not for extended discussion; this conversation has been moved to chat. Jul 21, 2016 at 15:14
• Here's the comment that had the last timing note since comments are supposed to be where this is posted. @TomWenseleers: Just redid the timings and the Julia and R ones indeed now are in exactly the same range, save for the graphics part that is! Great to hear also that the new version will be able to get this kind of performance even with my original code (with the const... fixed, which is the easy part)! Julia could become a great replacement also for Rcpp when callable for R - I gotta try it, stackoverflow.com/questions/9965747/linking-r-and-julia Jul 21, 2016 at 15:20

in your MATLAB code, some unnecessary steps can be skipped. Here's a modification of your doStep function that should improve performance quite a bit. Each modif is explained in the comments within the code.

doStep2.m function file:

function [ newtv ] = doStep2( tv, fitnessf, mutrate, stdev )
% Function to calculate new offspring trait values from parent population trait values tv
% given a density-dependent fitness function fitnessf,
% a mutation rate mutrate and standard deviation of mutational effects stdev

n = numel(tv) ; % numel is faster

numberoffspring = poissrnd(fitnessf(tv, n)) ;
newtv = repelem(tv, numberoffspring);

n = numel(newtv) ;

% logical vector for which of the offspring will endergo a mutation
mutants = rand(1,n) < mutrate ;

% faster way of generating normally distributed random numbers
% see Stewie's link in his answer
rnoise = randn(1, sum(mutants)) .* stdev ;

newtv(mutants) = newtv(mutants) + rnoise ;

% Check only the mutants for negative trait values
newtv(mutants) = max(newtv(mutants),0) ;

% the output is now newtv. This way, the argument tv is not altered
% within the function. Hence, MATLAB does not need to overwrite it.

end


On my laptop (i7 2.6Ghz, Windows 10, MATLAB 2015a) and with ngens = 100000 and psize = 1000, the timings are:

• With doStep: 81.004 seconds
• With doStep2: 54.679 seconds

Still not great compared to R but it's better.

Edit

Here's a more detailed explanation. In the new version doStep2, the line mutants = rand(1,n) < mutrate is basically what the binornd function does internally (type edit binornd to check). However, binornd throws out this vector and only returns its sum. Here we keep this logical vector to have the position of each mutant. This then allows avoiding the datasample function, which unnecessarily generates additional random numbers to select the position of the mutants. To summarize, the mutants' position is computed only once, not twice.

• Could you explain why your version of the function is faster? That will improve your answer. Aug 23, 2016 at 14:27
• Sure @pacmaninbw. I edited the response to add a more detailed explanation. Aug 23, 2016 at 14:50
• Very nice answer! ^ from me =) Sep 14, 2016 at 5:38
• Also, please have a look at my answer. I see that you use repelem and have R2015a. The code snippet I have in my answer should improve performance, but I'm not sure about hos much. Not necessarily much for the vector sizes here, but it's worth a try. :) Sep 14, 2016 at 5:48

Regarding the rpois function, I think the preferred way of simulating a single value from a Poisson distribution with mean λ is

using Distributions
rand(Poisson(λ))


To fill out, a vector of Ints, with Poisson simulated values from a vector of means, p, you can use the mutating version of map. In v0.4 it is best to write a helper function

r1(λ) = rand(Poisson(λ))


then call it as map!(r1, out, p). In v0.5 you can just as easily use an anonymous function in the map map!(λ -> rand(Poisson(λ)), out, p).

• We tested rpois using both StatsFuns.RFunctions.poisrand(p[i]) and rand(Poisson(p[i])). The first ends up being a little faster since the second is actually syntactic sugar for the first as defined here. Jul 22, 2016 at 18:01
• But it won't always be defined that way. The objective is to replace the calls to the C functions in Rmath library with native Julia code. Jul 22, 2016 at 19:20
• Yes, and I'll update it soon depending on the answer given here. Likely it will go back to rand(Poisson(lambda)), but the method this will be different. Jul 22, 2016 at 20:54

Firstly, practical comments regarding Julia implementation. The line:

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


consumes a lot of time because fitnessfunc is a user defined function about which not enough is known during compilation. This can be fixed by adding a type assertion to the result thus:

numberoffspring =
[rand(Poisson(theta)) for theta in fitnessfunc(tv, n, k)::Vector{Float64}]


Chris Rackauckas gives many details and especially a Package with an implementation of R's rep function. This useful function will eventually find its way into a natural Julia form. Meanwhile, it might be simpler to add it to the source thus:

function repeach{T}(v::Vector{T},each::Vector{Int})
length(v)==length(each) || throw(ArgumentError("values and each vectors must have same size"))
n = length(each)
s = sum(each)
r = Array(T,s)
iv = 1
for i = 1:n
@inbounds val = v[i]
@inbounds for j=1:each[i]
@inbounds r[iv] = val
iv += 1
end
end
return r
end


This is slightly ugly, but as Chris showed, this will be hidden normally. A one-liner version of this function would be

repeach2(v,each) = vcat([[v[i] for j=1:each[i]] for i=1:length(v)]...)


But it is less efficient. The relevant line in the code which uses repeach can be:

newtv = repeach(tv, numberoffspring) # offspring are copies of their parents


(Incidentally, this line as written in the post had a typo for me using Chris' Package).

To get a feel for the running time change, the poster should run the modified version.

Lastly, these types of cross language benchmark questions are vague. The initial implementations are usually in favor of the initial expertise of writers. Optimizations usually make a big difference, but it is unclear what is expected to be natural to a typical user of a language. Perhaps a language benchmark should consider the performance of a newbie user? When experts code a problem, many times languages converge on a single language agnostic bottleneck (for example Poisson RNG in this case) and the code converges into a too-complex for novice version. To summarize, these comparisons are problematic, but still useful and educational.

P.S. Considering the responses to this question and the prominence of Julia in them, if you would rate a language by community enthusiasm and future, Julia scores points.

• Many thanks for this - will try it and update you on the timings. In this case I was indeed most after how each language would perform in the hands of beginners (a course with biology students), but using typical paradigms you would use in each language (e.g. using vectorised code when possible etc). I am quite enthusiastic too about Julia, but feel it may still take one or two years before it's really ready from deployment in the classroom. Jul 21, 2016 at 20:35
• @DanGetz "(Incidentally, this line as written in the post had a typo for me using Chris' Package)" are you saying there's a typo still on the repo? If you run Pkg.update() do you still have the error? I don't have unit testing for any of this yet so hopefully the repo code works on other systems like it does on my test system! Jul 22, 2016 at 0:36
• "in the post had a typo". The line is R.rep(tv, each = numberoffspring) which explicitly names the each parameter although it is an optional positional parameter in the repo. Perhaps the repo has changed since then. BTW using 0.5. Jul 22, 2016 at 14:20
• Yes there was a change to the repo. Dispatch is weird on keyword arguments, and I didn't want to have to handle that right now. Jul 22, 2016 at 20:55

MATLAB

I don't have the Statistical Toolbox, so unfortunately I can't do any benchmarking to test this, but I believe the following points will increase the performance of your implementation. It will still be far slower than R though.

My guess is that the functions that are part of the toolbox are the slow ones: normrnd, poissrnd (and possibly datasample) takes up most of the time.

To improve the performance of normrnd you might benefit from doing something along the lines of this answer on SO.

To improve the performance of poissrnd you might benefit from either creating a custom function from scratch, or call the function inside a loop with few elements in each iteration (yup, counter-intuitive). The function I linked to is old, so Mathworks have probably improved their own poissrnd function by now, but it's worth a try. I've seen the loop-approach used to improve performance a few places but I can't find any reference now.

I don't know which version of MATLAB you're running. If it's R2015a then repelem is not the fastest way to create the newtv-vector for a large number of elements. The following function can process twice as many elements per second if the data size becomes larger than 10000 elements.

function out = rle_cumsum_diff(vals,runlens)
clens = cumsum(runlens);
idx(clens(end))=0;
idx([1 clens(1:end-1)+1]) = diff([0 vals]);
out = cumsum(idx);
return;


Your evolve function can probably be improved by pre-allocating memory to the out-cell array (use NaN). Since you don't know the number of elements prior to the doStep-function call, you must overshoot a bit, and remove NaN-elements after the function call.

This is definitely not the bottleneck in your code, but for reference:

numel is faster than length

Use timeit when timing functions in Matlab, it's more accurate than tic/toc.

• Many thanks for this good advice - I'll try it sooner or later! Jul 26, 2016 at 11:46