R Triangular Mesh

Question

I'm in the process of doing some development for the ggtern package, and I am trying to produce an efficient algorithm to do with ternary heatmaps. Specifically, I am using the following post (Ternary Heatmap) as a startingpoint.

Consider the function below, which is based off (part of) the above link:

# Produce Points for Triangular Mesh
triMesh = function( n = 1){
  n    = max(as.integer(n[1]),1)
  p    = data.frame()
  cnt  = 0
  inc  = 1.0/n
  stp  = seq(0,1,length.out = n + 1 )

  for (ix in seq_along(stp)){
    z <- stp[ix]
    y <- 1 - z
    x <- 0
    while ( y >= 0 ) {
      p   <- rbind(p, c(cnt, x, y, z))
      y   <- y - inc #Step x down
      x   <- x + inc #Step y up
      cnt <- cnt + 1          #Increment Count
    }
  }
  colnames(p) = c("IDPoint","x","y","z")
  p = round(p[with(p,order(y,x,-z)),],6)
  rownames(p) = 1:nrow(p) - 1
  p
}

AND here is my version, which is syntactically much more concise:

# Produce Points for Triangular Mesh
triMesh2 = function( n = 1 ){
  n   = as.integer(max(n[1],1))

  #Nested plyr calls
  result = ldply(0:n,function(y){ ##OUTER
    ldply(0:(n-y),function(x){    ##INNER
      data.frame(x,y,z = n -x -y) ##DIFF
    })
  })

  result        = data.frame( 1:nrow(result)-1,result/n)
  names(result) = c('IDPoint','x','y','z')
  result
}

Now, using microbenchmark the first algorithm finishes several times faster:

> microbenchmark(triMesh(10))
Unit: milliseconds
        expr      min      lq     mean   median       uq      max neval
 triMesh(10) 6.447525 6.91798 8.432698 7.334905 8.727805 23.37242   100

> microbenchmark(triMesh2(10))
Unit: milliseconds
         expr      min       lq     mean   median       uq     max neval
 triMesh2(10) 27.26659 29.34891 32.50808 31.43524 34.92925 51.8585   100
>

I was wondering if anyone can improve the performance of the second algorithm to something in the vicinity of the first one (or better)...

Cheers

minem · Answer 1 · 2018-01-18T15:06:38.303

Usually it is mush faster to simply use vectors :

triMesh3 <- function(n = 1){
  n <- as.integer(max(n[1], 1))
  result <- lapply(0:n, function(y){ 
    l <- lapply(0:(n - y), function(x){    
      c(x = x, y = y, z = n - x - y) 
    })
    Reduce(rbind, l)
  })
  result <- Reduce(rbind, result)
  row.names(result) <- NULL
  result <- cbind(1:nrow(result) - 1, result/n)
  result <- as.data.frame(result)
  names(result) <- c('IDPoint', 'x', 'y', 'z')
  result
}

all.equal(triMesh3(12), triMesh2(12))
# [1] TRUE
microbenchmark::microbenchmark(triMesh3(10),
                               triMesh2(10), 
                               triMesh(10), times = 100, unit = "relative")
# Unit: relative
#         expr      min       lq     mean   median       uq      max neval cld
# triMesh3(10)  1.00000  1.00000  1.00000  1.00000  1.00000  1.00000   100 a  
# triMesh2(10) 92.16829 89.07131 86.66111 88.32173 85.68915 63.29785   100   c
#  triMesh(10) 30.60108 29.70537 29.61635 29.83430 30.11924 32.40393   100  b

@CPak You beat me a little bit, I also wanted to update my answer:

triMesh4_minem <- function(n = 1){
  n <- as.integer(max(n[1], 1))
  y1 <- 0:n
  ys <- n - y1 + 1
  y <- sapply(1:(n + 1), function(x) y1[1:ys[x]])
  y <- unlist(y)
  x <- rep(y1, times = ys)
  result2 <- cbind(1:(length(x)) - 1, y/n, x/n, (n - y - x)/n)
  result <- as.data.frame.matrix(result2)
  names(result) <- c('IDPoint', 'x', 'y', 'z')
  result
}
all.equal(triMesh4_minem(2), triMesh4_cpack(2))
# [1] TRUE

microbenchmark::microbenchmark(triMesh4_minem(1e4),
                               triMesh4_cpack(1e4),
                               times = 10, unit = "relative")
# Unit: relative
#                 expr      min       lq     mean   median       uq      max neval cld
# triMesh4_minem(10000) 2.659507 2.572209 2.121967 1.965973 1.906203 1.905907    10   b
# triMesh4_cpack(10000) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    10  a

thanks very much, even though you got to a similar result in the end, due to the scaling issue with N, I have given it to CPak. — Nicholas Hamilton, Jan 18 '18 at 18:41

score 1 · Accepted Answer · answered Jan 18 '18 at 14:51

I wanted to provide an alternative way of generating your data. The issue with Reduce(...) is that it generally does not scale well with increasing N

triMesh4 <- function(n=1) {
    n <- as.integer(max(n[1], 1))
    temp <- seq(0, n, 1)
    df <- data.frame(
                x = unlist(sapply((n+1):1, function(i) temp[1:i])),
                y = rep(0:n, (n+1):1)
            )
    df$z <- n - df$x - df$y
    df <- cbind(0:(nrow(df)-1), df / n)
    names(df) <- c('IDPoint', 'x', 'y', 'z')
    return(df)
}

all.equal(triMesh3(12), triMesh4(12))
# [1] TRUE

library(microbenchmark)
N <- c(12, 16, 100)
lapply(N, function(i) microbenchmark(triMesh3(i), triMesh4(i), times=10L, unit="relative"))
# [[1]]
# Unit: relative
        # expr      min       lq     mean   median       uq      max neval
 # triMesh3(i) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    10
 # triMesh4(i) 1.484984 1.472767 1.466758 1.474142 1.470987 1.392629    10

# [[2]]
# Unit: relative
        # expr      min       lq     mean   median       uq       max neval
 # triMesh3(i) 1.075225 1.081014 1.017441 1.024083 1.015504 0.8398393    10
 # triMesh4(i) 1.000000 1.000000 1.000000 1.000000 1.000000 1.0000000    10

# [[3]]
# Unit: relative
        # expr      min       lq     mean  median       uq      max neval
 # triMesh3(i) 23.67992 23.33367 22.79632 23.2149 21.89245 21.32084    10
 # triMesh4(i)  1.00000  1.00000  1.00000  1.0000  1.00000  1.00000    10

R Triangular Mesh

2 Answers2