Skip to content

Instantly share code, notes, and snippets.

@pao

pao/PKbench.jl

Created March 29, 2012 20:00
Show Gist options
  • Save pao/2243086 to your computer and use it in GitHub Desktop.
Save pao/2243086 to your computer and use it in GitHub Desktop.
Download ZIP
Benchmarking some Julia functions derived from R functions
Raw
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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment