54 of 66 deleted 16 characters in body

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

enter image description here

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 = R.rpois(n,fitnessfunc(tv, n))
    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 = 100000  # 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)

@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")
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, xlabel = "Reproductive rate R");
plot(p1, p2)

enter image description here

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] 

enter image description here

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)

enter image description here

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 : 4.8 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 : 48 s
Matlab 2015a : 328 s
Mathematica 10 : 966 s

From this, R and Julia would appear to be the clear winners, and Julia would almost certainly win over R for applications with a lot of for...to loops etc, where R usually tends to perform quite poorly. 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.

EDIT: updated my Julia timings above based on the suggested optimisations below (many thanks for that!). I'll still wait a bit to award my +50 bounty to see what anyone else might come up with!