Scientific programming

SciProg2020.md

What makes a great scientific computing language?

• Expressivity: should be possible to express mathematical & programming concepts concisely & elegantly. This is important to reduce the mental overhead required for writing down and reading models and algorithms in the form of code.
• Speed: should be possible to convert code into software/program that gets the most out of modern computer hardware
  • Functional Programming
  • Parallel computing
• Extensibility: should have a powerful, unambiguous, easy to use package management system built right into the core of the language
• Visualisation: a large and very important part of scientific programming is generating side-effects: writing reports, generating plots, creating interactive documents

Why R?

• Statistical programming language
• Open source
• Good package management (CRAN, BioconductoR)
• Supports functional programming
• Some R code compact, easy to read and maintain

Why not R?

• Not well designed or easy to learn
• Vectorisation: sacrificed performance of for loops, not always implemented (sapply, lapply etc.)
• Package management excellent, but badly split (CRAN vs. BioconductoR)
• Conflicting and contrasting coding styles (base R vs. tidyverse vs. RCpp)
• Uses multiple CPUs in completely different ways on Linux vs. Windows, two competing packages (parallel & multicore)

Why Python?

• General purpose computing language
• Open source
• Exceptionally clear code, easy to learn
• Excellent package management (PIP)
• Partially supports functional programming

Why not Python?

• Scientific computing capability solely from other packages (numpy, pandas, scipy, cython): different languages to learn
• Not well designed for computing with multiple CPUs

Why Julia?

• Scientific programming language
• Open source
• Excellent package management, based on gitHub
• Supports functional programming
• Designed to take advantage of computers with multiple CPUs

Why not Julia?

• Indeed...

Example code

R: Plotting a histogram

synthetic_data = rnorm(10000,0,1)
hist(synthetic_data)

julia: Plotting a histogram

using Plots, Distributions
synthetic_data = rand(Normal(0,1),10000);
histogram(synthetic_data)

R: Simulate rolling a dice a large number of times

dicevals = 1:6
Nrolls = 10000000
randrolls = sample(dicevals, size=Nrolls, replace=TRUE)

R: Simulate dice rolls using 12 CPUs

library(parallel)
cl = makeCluster(12)
clusterExport(cl, varlist=c("dicevals"))
randrolls = parSapply(cl, 1:Nrolls, function(i) sample(dicevals, size=1))
stopCluster(cl)

julia: Simulate rolling a dice a large number of times

dicevals = 1:6;
Nrolls = 10000000;
randrolls = rand(dicevals, Nrolls);

julia: Simulate dice rolls using 12 CPUs

randrolls = zeros(Int8, Nrolls);
@inbounds Threads.@threads for i in 1:Nrolls
    randrolls[i] = rand(dicevals)
end