R - 加速嵌套循环(矢量化？),调整不同大小的序列

Question

对于具有3列(X位置,Y位置和某些值VAL)的数据集,我想对包含在某些XxY间隔/区间中的所有VAL执行一些操作(例如,平均值)(即,我想对我的空间进行网格化) ).

我最初写的幼稚下面的函数来这样做(myT是通过数据集,xbounds并且ybounds是连续的间隔中断点(箱)的向量):

calcPerBin1 <- function(myT, xbounds, ybounds) {
  newT <- data.frame(matrix(0, nrow=(length(xbounds)-1)*(length(ybounds)-1), ncol=3))
  names(newT) <- c("X","Y","MEAN")
  line <- 1
  for (i in 1:(length(xbounds)-1)) {
    for (j in 1:(length(ybounds)-1)) {
      myTsubset <- myT[myT$X >= xbounds[i] & myT$X < xbounds[i+1] &
                       myT$Y >= ybounds[j] & myT$Y < ybounds[j+1], ]
      newT$MEAN[line] <- mean(myTsubset$VAL)
      newT$X[line] <- mean(c(xbounds[i], xbounds[i+1]))
      newT$Y[line] <- mean(c(ybounds[j], ybounds[j+1]))
      line <- line+1
    }
  }
  return(newT)
}

SHORTCUT问题:如何改进上面的代码？(以下是我的第一次尝试 - 如果太长时间可以跳过!)

---

双for循环当然非常不理想,它的执行时间很糟糕(无法用我的真实数据集).我因此尝试了以下代码(即,如果我没有错误,内部循环被矢量化):

calcPerBin2 <- function(myT,xbounds, ybounds) {
  newT <- data.frame(matrix(0, nrow=(length(xbounds)-1)*(length(ybounds)-1), ncol=3))
  names(newT) <- c("X","Y","MEAN")
  xboundsmean <- vector() ; yboundsmean <- vector()
  for (i in 1:(length(xbounds)-1)) {
    xboundsmean <- c(xboundsmean, mean(c(xbounds[i],xbounds[i+1])))}
  for (i in 1:(length(ybounds)-1)) {
    yboundsmean <- c(yboundsmean, mean(c(ybounds[i],ybounds[i+1])))}
  xyvals <- expand.grid(xmid=xboundsmean, ymid=yboundsmean)
  xyvals$xmin <- xyvals$xmid-binsize/2
  xyvals$xmax <- xyvals$xmid+binsize/2
  xyvals$ymin <- xyvals$ymid-binsize/2
  xyvals$ymax <- xyvals$ymid+binsize/2
  res <- vector()
  for (i in 1:dim(xyvals)[1]) {
    cond <- (myT$X >= xyvals$xmin[i] & myT$X < xyvals$xmax[i] &
               myT$Y >= xyvals$ymin[i] & myT$Y < xyvals$ymax[i])
    res <- c(res, mean(myT$VAL[cond]))
  }
  newT$MEAN <- res
  newT$X <- xyvals[,1]
  newT$Y <- xyvals[,2]
  return(newT)
}

这非常难看,所以我尝试了以下变体:

calcPerBin2.2 <- function(myT,xbounds, ybounds, sizeofbin) {
  newT <- data.frame(matrix(0, nrow=(length(xbounds)-1)*(length(ybounds)-1), ncol=3))
  names(newT) <- c("X","Y","MEAN")
  xcut <- cut(myT$X, breaks=xbounds)
  ycut <- cut(myT$Y, breaks=ybounds)
  xycut <- expand.grid(XCUT=levels(xcut), YCUT=levels(ycut))
  xylowers <- cbind(xlower = as.numeric(sub("\\((.+),.*", "\\1", xycut$XCUT) ),
                    ylower = as.numeric(sub("\\((.+),.*", "\\1", xycut$YCUT) ))
  res <- vector()
  for (i in 1:dim(xycut)[1]) {
    cond <- (xcut==xycut$XCUT[i] & ycut==xycut$YCUT[i])
    res <- c(res, mean(myT$VAL[cond]))
  }
  newT$MEAN <- res
  newT$X <- xylowers[,1]+sizeofbin/2
  newT$Y <- xylowers[,2]+sizeofbin/2
  return(newT)
}

我可以运行它,例如:

# Control parameters
xmax <- 500 
ymax <- 1000 
N <- 100000
binsize <- 50

xbins <- seq(0,xmax,binsize)
ybins <- seq(0,ymax,binsize) # xbins and ybins do NOT have the same size

# Generate dummy data
xcoords <- runif(N, 1, xmax) 
ycoords <- runif(N, 1, ymax) 
vals <- xcoords+ycoords**2
data <- data.frame(cbind(X=xcoords, Y=ycoords, VAL=vals))

# Run
system.time(test1 <- calcPerBin1(data, xbins, ybins))
system.time(test2 <- calcPerBin2(data, xbins, ybins))
system.time(test2.2 <- calcPerBin2.2(data, xbins, ybins, binsize))

产生了轻微的改善calcPerBin2,但calcPerBin2.2甚至更糟calcPerBin1(并且,是的,所有代码都很难看).我的问题是,我不太清楚如何替换(矢量化？)剩余的循环calcPerBin2.例如,如何在矢量形式(它们不具有相同的大小)而不是我使用的索引形式myT$X的基础上有效地写入条件？任何改进上述代码的建议都是受欢迎的 - 谢谢.xyvals$xmincalcPerBin2

Answer 1

您可以在三行中完成大部分操作(使用zoofor rollmean):

library(zoo) # load the package
data$X <- cut(data$X, xbins, labels = rollmean(xbins, 2))
data$Y <- cut(data$Y, ybins, labels = rollmean(ybins, 2))
res <- aggregate(VAL ~ X + Y, data, mean)

检查结果:

# order it the same way as in test1, then show the first lines
head(res[order(res$X, res$Y),]) 
#    X   Y        VAL
#1  25  25   900.8305
#11 25  75  5957.4972
#21 25 125 15680.8103
#31 25 175 30877.6696
#41 25 225 50688.4860
#51 25 275 75961.8558

将其与原始函数的结果进行比较:

test1 <- calcPerBin1(data, xbins, ybins)
head(test1)
#   X   Y       MEAN
#1 25  25   900.8305
#2 25  75  5957.4972
#3 25 125 15680.8103
#4 25 175 30877.6696
#5 25 225 50688.4860
#6 25 275 75961.8558

基准测试:

fastbin <- function(data, xbins, ybins){
  data$X <- cut(data$X, xbins, labels = rollmean(xbins, 2))
  data$Y <- cut(data$Y, ybins, labels = rollmean(ybins, 2))
  aggregate(VAL ~ X + Y, data, mean)
}

library(dplyr)   # for faster aggregation
fastbin.dplyr <- function(data, xbins, ybins){
  data %>%
    mutate(X = cut(X, xbins, labels = rollmean(xbins, 2)),
           Y = cut(Y, ybins, labels = rollmean(ybins, 2))) %>%
    group_by(X, Y) %>% 
    summarise(Val = mean(VAL))
}

system.time(test1 <- calcPerBin1(data, xbins, ybins))
       User      System     elapsed 
       3.47        0.12        3.59 

system.time(res.fastbin <- fastbin(data, xbins, ybins))
       User      System     elapsed 
       1.01        0.02        1.05 

system.time(res.fastbin.dplyr <- fastbin.dplyr(data, xbins, ybins))
       User      System     elapsed 
       0.06        0.00        0.06