加快茱莉亚写得不好的R例子
Julia比较性能和R的例子看起来特别复杂 。 https://github.com/JuliaLang/julia/blob/master/test/perf/perf.R
什么是最快的性能，你可以从下面的两个algorithm中得出（最好是解释一下你改变了什么使它更像R）？
## mandel mandel = function(z) { c = z maxiter = 80 for (n in 1:maxiter) { if (Mod(z) > 2) return(n-1) z = z^2+c } return(maxiter) } mandelperf = function() { re = seq(-2,0.5,.1) im = seq(-1,1,.1) M = matrix(0.0,nrow=length(re),ncol=length(im)) count = 1 for (r in re) { for (i in im) { M[count] = mandel(complex(real=r,imag=i)) count = count + 1 } } return(M) } assert(sum(mandelperf()) == 14791) ## quicksort ## qsort_kernel = function(a, lo, hi) { i = lo j = hi while (i < hi) { pivot = a[floor((lo+hi)/2)] while (i <= j) { while (a[i] < pivot) i = i + 1 while (a[j] > pivot) j = j - 1 if (i <= j) { t = a[i] a[i] = a[j] a[j] = t } i = i + 1; j = j - 1; } if (lo < j) qsort_kernel(a, lo, j) lo = i j = hi } return(a) } qsort = function(a) { return(qsort_kernel(a, 1, length(a))) } sortperf = function(n) { v = runif(n) return(qsort(v)) } sortperf(5000) 嗯，在Mandelbrot的例子中，matrixM的尺寸是转置的
M = matrix(0.0,nrow=length(im), ncol=length(re))
因为它是通过在内部循环中增加count来填充的（ im连续值）。 我的实现在mandelperf.1创build了一个复数的向量，并使用索引和子集对所有元素进行操作，以跟踪向量中哪些元素尚未满足条件Mod(z) <= 2
mandel.1 = function(z, maxiter=80L) { c <- z result <- integer(length(z)) i <- seq_along(z) n <- 0L while (n < maxiter && length(z)) { j <- Mod(z) <= 2 if (!all(j)) { result[i[!j]] <- n i <- i[j] z <- z[j] c <- c[j] } z <- z^2 + c n <- n + 1L } result[i] <- maxiter result } mandelperf.1 = function() { re = seq(-2,0.5,.1) im = seq(-1,1,.1) mandel.1(complex(real=rep(re, each=length(im)), imaginary=im)) }
为13倍加速（结果是相等的，但不相同，因为原来的返回数字而不是整数值）。
> library(rbenchmark) > benchmark(mandelperf(), mandelperf.1(), + columns=c("test", "elapsed", "relative"), + order="relative") test elapsed relative 2 mandelperf.1() 0.412 1.00000 1 mandelperf() 5.705 13.84709 > all.equal(sum(mandelperf()), sum(mandelperf.1())) [1] TRUE
快速sorting的例子实际上没有sorting
> set.seed(123L); qsort(sample(5)) [1] 2 4 1 3 5
但是我的主要加速是围绕主轴vector化分区
qsort_kernel.1 = function(a) { if (length(a) < 2L) return(a) pivot <- a[floor(length(a) / 2)] c(qsort_kernel.1(a[a < pivot]), a[a == pivot], qsort_kernel.1(a[a > pivot])) } qsort.1 = function(a) { qsort_kernel.1(a) } sortperf.1 = function(n) { v = runif(n) return(qsort.1(v)) }
对于7倍的加速（相对于未修正的原始）
> benchmark(sortperf(5000), sortperf.1(5000), + columns=c("test", "elapsed", "relative"), + order="relative") test elapsed relative 2 sortperf.1(5000) 6.60 1.000000 1 sortperf(5000) 47.73 7.231818
由于在原始比较中，Julia比mandel的R快30倍，而quicksort的快500倍，上面的实现仍然没有真正的竞争力。
这个问题的关键词是“algorithm”：
什么是最快的性能，你可以从下面的两个algorithm中得出（最好是解释一下你改变了什么使它更像R）？
就像“你能在R里面使用这些algorithm有多快？ 这里所讨论的algorithm是标准的Mandelbrot复杂循环迭代algorithm和标准的recursion快速sorting内核。
计算这些基准testing中出现的问题的答案当然有更快的方法 – 但不使用相同的algorithm。 您可以避免recursion，避免迭代，并避免R不擅长的任何其他。 但是，你不再比较相同的algorithm。
如果你真的想要用R来计算Mandelbrot集合或者对数字进行sorting，那么这不是你写代码的方式。 您可以尽可能多地对其进行vector化，从而将所有工作都推送到预定义的C内核中，或者只是编写一个自定义的C扩展并在那里进行计算。 无论哪种方式，结论是R不够快，以获得真正的好performance – 你需要C做大部分的工作，以获得良好的performance。
这正是这些基准testing的要点：在Julia，你永远不必依靠C代码来获得良好的性能。 你可以写纯粹的朱莉娅想做的事情，并且会有很好的performance。 如果一个迭代标量循环algorithm是做自己想做的最自然的方法，那么就这样做。 如果recursion是解决问题的最自然的方法，那也没关系。 在任何时候，你都不得不依靠C来获得性能 – 无论是通过非自然的vector化还是编写自定义的C扩展。 当然，你可以写vector化的代码，因为它通常是线性代数; 你可以打电话给C，如果你已经有一些图书馆做你想做的事情。 但是你不需要。
我们希望在不同的语言中使用相同的algorithm进行最公平的比较：
- 如果有人在R中使用相同的algorithm有更快的版本，请提交补丁！
- 我相信在julia网站上的R基准testing已经是字节编译的，但是如果我做错了，比较对R来说是不公平的，请告诉我，我会修复它并更新基准。
