55 of 66 Rollback to Revision 53

# Benchmarks of scientific programming languages (R, Julia, Mathematica, Matlab) for theoretical modelling in biology course

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 re optimising code performance (say in R using Rcpp, data.table or dplyr), 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, although an elegant C, C++, JAVA or JVM OO-based programming language implementation would be cool too. Any help with this would be much appreciated, as it would help me to decide on the best platform for my classes!

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
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)[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(tv, fitnessfunc, mutrate, stdev, ngens)
Results = Array(Array, ngens)
Results[1] = 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)[2]/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)? Other elegant solutions, possibly in other languages, that would scale better for larger problems (ngens or pop sizes) are also welcome! I'll award my +50 bounty to whoever can come up with a solution that strikes the best balance between compactness/elegance/readability and performance and scalability. 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.