dirkschumacher/armacmp_julia.md

## armacmp_julia.md

      
    Raw
  

              armacmp_julia.md
            
          
    library(armacmp)

# some of julia's microbenchmarks translated to C++
# https://github.com/JuliaLang/Microbenchmarks/blob/master/perf.R

fib_cpp <- armacmp(function(n = type_scalar_int()) {
  fib_rec <-  function(nr = type_scalar_int()) {
    if (nr < 2) {
      return(nr, type = type_scalar_int())
    } else {
      return(fib_rec(nr-1) + fib_rec(nr-2), type = type_scalar_int())
    }
  }
  return(fib_rec(n), type = type_scalar_int())
})
stopifnot(fib_cpp(20) == 6765)
fib_r <- function(n) {
  fib_rec <-  function(nr) {
    if (nr < 2) {
      return(nr)
    } else {
      return(fib_rec(nr-1) + fib_rec(nr-2))
    }
  }
  return(fib_rec(n))
}
stopifnot(fib_r(20) == 6765)

microbenchmark::microbenchmark(
  fib_cpp(20),
  fib_r(20)
)
#> Unit: microseconds
#>         expr      min         lq      mean     median        uq      max
#>  fib_cpp(20)   66.984    77.4265   227.237    81.5845   116.228 12760.85
#>    fib_r(20) 7912.313 10791.2160 16243.916 14095.6780 17064.966 69667.72
#>  neval
#>    100
#>    100

pisum_cpp <- armacmp(function() {
  t <- 0.0
  for (j in seq_len(500L)) {
    t <- 0.0
    for (k in seq_len(10000L)) {
      t <- t + 1.0/(k*k)
    }
  }
  return(t, type = type_scalar_numeric())
})
stopifnot(abs(pisum_cpp()-1.644834071848065) < 1e-12)

pisum_r <- function() {
  t <- 0.0
  for (j in seq_len(500L)) {
    t <- 0.0
    for (k in seq_len(10000L)) {
      t <- t + 1.0/(k*k)
    }
  }
  return(t)
}
stopifnot(abs(pisum_r()-1.644834071848065) < 1e-12)

microbenchmark::microbenchmark(
  pisum_cpp(),
  pisum_r()
)
#> Unit: milliseconds
#>         expr        min         lq      mean    median         uq
#>  pisum_cpp()   7.432039   7.666541  10.21577   8.07064   9.061706
#>    pisum_r() 266.758123 329.065431 387.83702 357.52156 394.785481
#>         max neval
#>    49.91694   100
#>  1201.24953   100

# julia microbenchmark qsort
qsort_r = function(a) {
  qsort_kernel = function(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(lo, j)
      lo = i
      j = hi
    }
  }
  qsort_kernel(1, length(a))
  return(a)
}

qsort_cpp <- armacmp(function(vec = type_colvec()) {
  vec_copy <- vec
  qsort_kernel <- function(lo = type_scalar_int(), hi = type_scalar_int()) {
    i <- lo
    j <- hi
    while (i < hi) {
      pivot <- vec_copy[floor((lo+hi)/2)]
      while (i <= j) {
        while (vec_copy[i] < pivot) i <- i + 1
        while (vec_copy[j] > pivot) j <- j - 1
        if (i <= j) {
          t <- vec_copy[i]
          vec_copy[i] <- vec_copy[j]
          vec_copy[j] <- t
          i <- i + 1
          j <- j - 1
        }
      }
      if (lo < j) qsort_kernel(lo, j)
      lo <- i
      j <- hi
    }
  }
  qsort_kernel(1, length(vec_copy))
  return(vec_copy, type = type_colvec())
})

x <- runif(10000)
all.equal(
  as.numeric(qsort_cpp(x)),
  qsort_r(x)
)
#> [1] TRUE
microbenchmark::microbenchmark(
  as.numeric(qsort_cpp(x)),
  qsort_r(x)
)
#> Unit: microseconds
#>                      expr       min         lq       mean     median
#>  as.numeric(qsort_cpp(x))   778.583   794.0085   889.9126   825.2585
#>                qsort_r(x) 36767.078 39312.2270 46781.8513 40993.8765
#>         uq       max neval
#>    923.019  1790.045   100
#>  52073.695 81240.261   100
^{Created on 2019-08-08 by the reprex package (v0.3.0)}