R中的矢量化索引答案

【问题标题】：Vectorizing Indexing in RR中的矢量化索引
【发布时间】：2015-09-30 22:51:18
【问题描述】：

我正在尝试在大型数据框中建立索引。 sdata 框架有 300 万个观测值和 26 个变量（底部）。

setsize <- 6
eccent <- 150
ctrX <- 400
ctrY <- 300

xyrotate <- function(x,y,ctrX,ctrY,angle){
  distX <- x - ctrX;
  distY <- y - ctrY;
  radians <- angle * (pi/180);
  rotX <- ctrX + (distX*cos(radians)) - (distY*sin(radians));
  rotY <- ctrY + (distX*sin(radians)) + (distY*cos(radians));
  coordinates <- list("X" = rotX,"Y" = rotY)
  return(coordinates)
}

loc <- data.frame(x = numeric(setsize),
              y = numeric(setsize))
loc$x[1] <- ctrX
loc$y[1] <- ctrY - eccent
for(i in 2:setsize){
  coord <- xyrotate(loc$x[1], loc$y[1],ctrX,ctrY,(i-1)*(360/setsize)) 
  loc$x[i] <- coord$X
  loc$y[i] <- coord$Y
}
gazedist <- matrix(nrow=nrow(sdata), ncol = setsize)
for(d in 1:setsize){
  x <- sdata$RIGHT_GAZE_X-loc$x[d]
  y <- sdata$RIGHT_GAZE_Y-loc$y[d]
  gazedist[,d] <- sqrt(x^2+y^2)
}
sdata$gdist_T <- 0
sdata$gdist_T <- gazedist[ ,sdata$t_targLoc]

这里的最后一行导致 R 崩溃。有没有办法将 sdata$t_targLoc[i] 的值插入到gazedistance[i,d] 的 d 值中。 for 循环等效项是：

for(i in 1:length(gazedist)){
   sdata$gdist_T[i] <- gazedist[i,sdata$t_targLoc[i]]
}

但这在 R 中会很慢...

s数据结构

structure(list(RIGHT_GAZE_X = c(409.5, 409.6, 409.5, 409.4, 409.3, 
409.2, 409.1, 409, 408.9, 408.8), RIGHT_GAZE_Y = c(291.9, 291.5, 
290.9, 290.3, 290.3, 290.3, 289.8, 289.2, 288.7, 288.8), RECORDING_SESSION_LABEL = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "ET101", class = "factor"), 
    t_block = c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), t_trialNum = c(129L, 
    129L, 129L, 129L, 129L, 129L, 129L, 129L, 129L, 129L), t_subjNum = c(101L, 
    101L, 101L, 101L, 101L, 101L, 101L, 101L, 101L, 101L), t_colCond = structure(c(1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "green", class = "factor"), 
    t_targLoc = c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), t_targID = structure(c(1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "L", class = "factor"), 
    t_targShape = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L), .Label = "diamond", class = "factor"), t_singLoc = c(5L, 
    5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), t_singPres = structure(c(1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "P", class = "factor"), 
    t_singDist = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), t_singAngle = c(120L, 
    120L, 120L, 120L, 120L, 120L, 120L, 120L, 120L, 120L), t_targAngle = c(120L, 
    120L, 120L, 120L, 120L, 120L, 120L, 120L, 120L, 120L), t_RESP = structure(c(1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "L", class = "factor"), 
    t_ACC = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), t_RT = c(686.1062, 
    686.1062, 686.1062, 686.1062, 686.1062, 686.1062, 686.1062, 
    686.1062, 686.1062, 686.1062), TRIAL_START_TIME = c(1031031L, 
    1031031L, 1031031L, 1031031L, 1031031L, 1031031L, 1031031L, 
    1031031L, 1031031L, 1031031L), TIMESTAMP = c(1031030, 1031032, 
    1031034, 1031036, 1031038, 1031040, 1031042, 1031044, 1031046, 
    1031048), IP_START_TIME = c(1031031L, 1031031L, 1031031L, 
    1031031L, 1031031L, 1031031L, 1031031L, 1031031L, 1031031L, 
    1031031L), currtime = c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18
    ), currsamp = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), gdist_T = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0), gdist_S = c(0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0), gdist_NS = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)), .Names = c("RIGHT_GAZE_X", 
"RIGHT_GAZE_Y", "RECORDING_SESSION_LABEL", "t_block", "t_trialNum", 
"t_subjNum", "t_colCond", "t_targLoc", "t_targID", "t_targShape", 
"t_singLoc", "t_singPres", "t_singDist", "t_singAngle", "t_targAngle", 
"t_RESP", "t_ACC", "t_RT", "TRIAL_START_TIME", "TIMESTAMP", "IP_START_TIME", 
"currtime", "currsamp", "gdist_T", "gdist_S", "gdist_NS"), row.names = 53170:53179, class = "data.frame")

【问题讨论】：

您意识到您不必上传所有行，可能只上传前 10 行？
"这里的最后一行导致 R 崩溃。" -> 有什么错误信息？该行本身是可以的，完全等同于循环（所以要么工作要么失败）。

标签： r

【解决方案1】：

您似乎想要获取gazedist 的第i 行和sdist$$t_targLoc[i] 列。有一个内置的。使用：

sdata$gdist_T <- gazedist[cbind(1:nrow(gazedist),sdata$t_targLoc)]

这是一个例子：

m <- matrix(1:25,nc=5)
m
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    1    6   11   16   21
# [2,]    2    7   12   17   22
# [3,]    3    8   13   18   23
# [4,]    4    9   14   19   24
# [5,]    5   10   15   20   25
v <- c(1,3,5,2,4)
m[cbind(1:5,v)]
# [1]  1 12 23  9 20

顺便说一句，您的第一个 for 循环可以替换为：

loc <- as.data.frame(xyrotate(ctrX,ctrY-eccent, ctrX,ctrY,(1:(setsize-1))*360/setsize))
loc <- rbind(c(X=ctrX,Y=ctrY-eccent),lloc)

您的第二个 for 循环可以替换为，例如，

f <- function(x,y) {
  x <- sdata$RIGHT_GAZE_X-x
  y <- sdata$RIGHT_GAZE_Y-y
  sqrt(x^2+y^2)
}
ggazedist <- mapply(f, loc$x, loc$y)
identical(gazedist,ggazedist)
# [1] TRUE

【讨论】：

【解决方案2】：

是的，我认为这行得通。本质上，您想从矩阵中获取一个向量，其中行索引为1:nrow(sdat)，列索引为sdata$t_targLoc。这不是内置的（据我所知），但我们可以将矩阵转换为向量并获取正确的值。

gazedist_vals = as.vector(gazedist)
rows = 1:nrow(sdat)
cols = sdat$t_targLoc
indices = (cols - 1) * nrow(gazedist) + rows
sdata$gdist_T = gazedist_vals[indices]

我认为这会满足你的要求。

您的数据摘录不能很好地说明这一点，因为t_targLoc 始终为 3。这里有一个小例子：

x = matrix(c(5,2,65,8,4,2), nrow = 2)
x
#      [,1] [,2] [,3]
# [1,]    5   65    4
# [2,]    2    8    2
as.vector(x)
# [1]  5  2 65  8  4  2
rows = c(1, 1, 2)
cols = c(3,2,1)
inds = (cols - 1) * nrow(x) + rows
as.vector(x)[inds]
# [1]  4 65  2

一条评论：您的问题很详尽，但最少示例通常更可取。您给了我们 26 列数据，其中只需要几列。你给了我们计算距离的函数和代码，而你本来可以只给出一个距离矩阵（然后只需要一个数据列）。像我的x 矩阵和rows 和cols 向量这样的小例子可能就是您展示问题所需要的全部内容。

【讨论】：

【解决方案3】：

你可以，但我怀疑你会节省很多时间，除非你找到一种完全矢量化的方法。换句话说，您必须避免使用像 apply 或 sapply 这样的函数，它们都基于 C 中的 for 循环，因此不会比正常的 for 循环快得多，如果有的话。

someFunction <- function(x) ifelse(x %in% seq(0, 50000, 100), 1, 0)

# Here you have "vectorized" the indexing
system.time(sapply(1:nrow(diamonds), someFunction))
# 2.6 elapsed secs

## vs here where you're just using a for loop

system.time(
for(i in 1:nrow(diamonds)) {
  k[i] <- someFunction(i)
}
)
# 2.7 elapsed secs

【讨论】：