Benchmarking some Julia functions derived from R functions

PKbench.jl

# in http://dmbates.blogspot.com/2012/03/pk-models-in-r-and-in-julia.html I describe
# the motivation behind these functions.  The test data is extracted from an R data set
# Use semicolons between elements to ensure this is a vector
Time = [0; 0.25; 0.57; 1.12; 2.02; 3.82; 5.1; 7.03; 9.05; 12.12; 24.37; 0; 0.27; 0.52; 1; 1.92; 3.5; 5.02; 7.03; 9; 12; 24.3; 0; 0.27; 0.58; 1.02; 2.02; 3.62; 5.08; 7.07; 9; 12.15; 24.17; 0; 0.35; 0.6; 1.07; 2.13; 3.5; 5.02; 7.02; 9.02; 11.98; 24.65; 0; 0.3; 0.52; 1; 2.02; 3.5; 5.02; 7.02; 9.1; 12; 24.35; 0; 0.27; 0.58; 1.15; 2.03; 3.57; 5; 7; 9.22; 12.1; 23.85; 0; 0.25; 0.5; 1.02; 2.02; 3.48; 5; 6.98; 9; 12.05; 24.22; 0; 0.25; 0.52; 0.98; 2.02; 3.53; 5.05; 7.15; 9.07; 12.1; 24.12; 0; 0.3; 0.63; 1.05; 2.02; 3.53; 5.02; 7.17; 8.8; 11.6; 24.43; 0; 0.37; 0.77; 1.02; 2.05; 3.55; 5.05; 7.08; 9.38; 12.1; 23.7; 0; 0.25; 0.5; 0.98; 1.98; 3.6; 5.02; 7.03; 9.03; 12.12; 24.08; 0; 0.25; 0.5; 1; 2; 3.52; 5.07; 7.07; 9.03; 12.05; 24.15]
# Or use the columnization indexing trick, this also works in MATLAB
#Time = Time[:]

Dose = [4.02; 4.02; 4.02; 4.02; 4.02; 4.02; 4.02; 4.02; 4.02; 4.02; 4.02; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 4.4; 5.86; 5.86; 5.86; 5.86; 5.86; 5.86; 5.86; 5.86; 5.86; 5.86; 5.86; 4; 4; 4; 4; 4; 4; 4; 4; 4; 4; 4; 4.95; 4.95; 4.95; 4.95; 4.95; 4.95; 4.95; 4.95; 4.95; 4.95; 4.95; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 4.53; 3.1; 3.1; 3.1; 3.1; 3.1; 3.1; 3.1; 3.1; 3.1; 3.1; 3.1; 5.5; 5.5; 5.5; 5.5; 5.5; 5.5; 5.5; 5.5; 5.5; 5.5; 5.5; 4.92; 4.92; 4.92; 4.92; 4.92; 4.92; 4.92; 4.92; 4.92; 4.92; 4.92; 5.3; 5.3; 5.3; 5.3; 5.3; 5.3; 5.3; 5.3; 5.3; 5.3; 5.3]

# This is creating the vectors of parameter values.  In practice they are not constant values.

lV  = fill(-1., length(Dose))
lka = fill(0.6, length(Dose))
lCl = fill(-4., length(Dose))

mf0(dose, t, V, ka, k) = (dose ./ V) .* (ka ./ (ka - k)) .* (exp(-k .* t) - exp(-ka .* t))
mf1(dose, t, lV, lka, lCl) = mf0(dose, t, exp(lV), exp(lka), exp(lCl - lV))

function mfg(dose, t, lV, lka, lCl)
    V      = exp(lV)
    expr2  = dose ./ V
    ka     = exp(lka)
    Cl     = exp(lCl)
    k      = Cl ./ V
    expr6  = ka - k
    expr7  = ka ./ expr6
    expr8  = expr2 .* expr7
    expr11 = exp(-k .* t)
    expr14 = exp(-ka .* t)
    expr15 = expr11 - expr14
    expr18 = V .* V
    expr19 = Cl .* V ./ expr18
    expr24 = expr6 .* expr6
    expr8 .* expr15, hcat(expr8 .* (expr11 .* (expr19 .* t)) - (expr2 .* (ka .* expr19 ./ expr24) + dose .* V ./ expr18 .* expr7) .* expr15,
                          expr2 .* (expr7 - ka .* ka ./ expr24) .* expr15 + expr8 .* (expr14 .* (ka .* t)),
                          expr2 .* (ka .* k ./ expr24) .* expr15 - expr8 .* (expr11 .* (k .* t)))
end

# I run one call to the functions to ensure that the method I will use is compiled

resA = mf1(Dose, Time, lV, lka, lCl)
resB = mfg(Dose, Time, lV, lka, lCl)[1]
# Check that these are closeish
@assert all(resA-resB < eps(max(max(resA), max(resB))))

@time for i=1:1000 mf1(Dose, Time, lV, lka, lCl) end
@time for i=1:1000 mfg(Dose, Time, lV, lka, lCl) end

# To run as batch, use julia PKbench.jl; I actually didn't know about -L, which will be nicer than
# putting a load() in the argument to julia -e, which is what I've been doing for a lot of tests.