gpt4 book ai didi

r - 在 R 中找到矩阵的相邻元素

转载 作者:行者123 更新时间:2023-12-04 13:16:07 25 4
gpt4 key购买 nike

编辑:非常感谢以下用户的巨大贡献和 Gregor 的基准测试。

假设我有一个填充了这样的整数值的矩阵......

    mat <- matrix(1:100, 10, 10)

我可以像这样创建每个元素的 x、y 坐标列表...
    addresses <- expand.grid(x = 1:10, y = 1:10)

现在对于这些坐标中的每一个(即对于 mat 中的每个元素),我想找到相邻元素(包括对角线,这应该是 8 个邻居)。

我确定有一个简单的方法,有人可以帮忙吗?

到目前为止,我所尝试的是循环并为每个元素记录相邻元素,如下所示;
    neighbours <- list()
for(i in 1:dim(addresses)[1]){
x <- addresses$x[i]
y <- addresses$y[i]
neighbours[[i]] <- c(mat[y-1, x ],
mat[y-1, x+1],
mat[y , x+1],
mat[y+1, x+1],
mat[y+1, x ],
mat[y+1, x-1],
mat[y , x-1],
mat[y-1, x-1])
}

当它碰到矩阵的边缘时会遇到问题,特别是当索引大于矩阵的边缘时。

最佳答案

这是一个很好的例子。我做了 4x4,所以我们可以很容易地看到它,但它可以通过 n 进行调整。 .它也是完全矢量化的,所以应该有很好的速度。

n = 4
mat = matrix(1:n^2, nrow = n)
mat.pad = rbind(NA, cbind(NA, mat, NA), NA)

使用填充矩阵,邻居只是 n × n 子矩阵,四处移动。使用罗盘方向作为标签:
ind = 2:(n + 1) # row/column indices of the "middle"
neigh = rbind(N = as.vector(mat.pad[ind - 1, ind ]),
NE = as.vector(mat.pad[ind - 1, ind + 1]),
E = as.vector(mat.pad[ind , ind + 1]),
SE = as.vector(mat.pad[ind + 1, ind + 1]),
S = as.vector(mat.pad[ind + 1, ind ]),
SW = as.vector(mat.pad[ind + 1, ind - 1]),
W = as.vector(mat.pad[ind , ind - 1]),
NW = as.vector(mat.pad[ind - 1, ind - 1]))

mat
# [,1] [,2] [,3] [,4]
# [1,] 1 5 9 13
# [2,] 2 6 10 14
# [3,] 3 7 11 15
# [4,] 4 8 12 16

neigh[, 1:6]
# [,1] [,2] [,3] [,4] [,5] [,6]
# N NA 1 2 3 NA 5
# NE NA 5 6 7 NA 9
# E 5 6 7 8 9 10
# SE 6 7 8 NA 10 11
# S 2 3 4 NA 6 7
# SW NA NA NA NA 2 3
# W NA NA NA NA 1 2
# NW NA NA NA NA NA 1

所以你可以看到第一个元素 mat[1,1] ,从北开始顺时针方向,邻居是 neigh 的第一列.下一个元素是 mat[2,1] ,依此类推 mat 的列. (您还可以与 @mrip 的答案进行比较,并看到我们的列具有相同的元素,只是顺序不同。)

基准测试

小矩阵
mat = matrix(1:16, nrow = 4)
mbm(gregor(mat), mrip(mat), marat(mat), u20650(mat), times = 100)
# Unit: microseconds
# expr min lq mean median uq max neval cld
# gregor(mat) 25.054 30.0345 34.04585 31.9960 34.7130 61.879 100 a
# mrip(mat) 420.167 443.7120 482.44136 466.1995 483.4045 1820.121 100 c
# marat(mat) 746.462 784.0410 812.10347 808.1880 832.4870 911.570 100 d
# u20650(mat) 186.843 206.4620 220.07242 217.3285 230.7605 269.850 100 b

在更大的矩阵上,我不得不取出 user20650 的函数,因为它试图分配一个 232.8 Gb 的向量,并且在等待大约 10 分钟后我也取出了 Marat 的答案。
mat = matrix(1:500^2, nrow = 500)

mbm(gregor(mat), mrip(mat), times = 100)
# Unit: milliseconds
# expr min lq mean median uq max neval cld
# gregor(mat) 19.583951 21.127883 30.674130 21.656866 22.433661 127.2279 100 b
# mrip(mat) 2.213725 2.368421 8.957648 2.758102 2.958677 104.9983 100 a

所以看起来在任何时间很重要的情况下,@mr​​ip 的解决方案是迄今为止最快的。

使用的功能:
gregor = function(mat) {
n = nrow(mat)
mat.pad = rbind(NA, cbind(NA, mat, NA), NA)
ind = 2:(n + 1) # row/column indices of the "middle"
neigh = rbind(N = as.vector(mat.pad[ind - 1, ind ]),
NE = as.vector(mat.pad[ind - 1, ind + 1]),
E = as.vector(mat.pad[ind , ind + 1]),
SE = as.vector(mat.pad[ind + 1, ind + 1]),
S = as.vector(mat.pad[ind + 1, ind ]),
SW = as.vector(mat.pad[ind + 1, ind - 1]),
W = as.vector(mat.pad[ind , ind - 1]),
NW = as.vector(mat.pad[ind - 1, ind - 1]))
return(neigh)
}

mrip = function(mat) {
m2<-cbind(NA,rbind(NA,mat,NA),NA)
addresses <- expand.grid(x = 1:4, y = 1:4)
ret <- c()
for(i in 1:-1)
for(j in 1:-1)
if(i!=0 || j !=0)
ret <- rbind(ret,m2[addresses$x+i+1+nrow(m2)*(addresses$y+j)])
return(ret)
}

get.neighbors <- function(rw, z, mat) {
# Convert to absolute addresses
z2 <- t(z + unlist(rw))
# Choose those with indices within mat
b.good <- rowSums(z2 > 0)==2 & z2[,1] <= nrow(mat) & z2[,2] <= ncol(mat)
mat[z2[b.good,]]
}

marat = function(mat) {
n.row = n.col = nrow(mat)
addresses <- expand.grid(x = 1:n.row, y = 1:n.col)
# Relative addresses
z <- rbind(c(-1,0,1,-1,1,-1,0,1), c(-1,-1,-1,0,0,1,1,1))
apply(addresses, 1,
get.neighbors, z = z, mat = mat) # Returns a list with neighbors
}

u20650 = function(mat) {
w <- which(mat==mat, arr.ind=TRUE)
d <- as.matrix(dist(w, "maximum", diag=TRUE, upper=TRUE))
# extract neighbouring values for each element
# extract where max distance is one
a <- apply(d, 1, function(i) mat[i == 1] )
names(a) <- mat
return(a)
}

关于r - 在 R 中找到矩阵的相邻元素,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/29105175/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com