如何使 R foreach 循环高效
How to make R foreach loops efficient
我正在尝试在 R 中计算一个 300,000x300,000 的矩阵，我的代码工作得很好，但现在已经 运行ning 好几天了，我怎样才能让它更高效、更省时？
我的代码运行良好，但已经 运行ning 好几天了，附件是我正在使用的代码的一个子集，ID 扩展到 300,000；我怎样才能使代码 运行 在几分钟内更快，因为它已经 运行 好几天了。
fam <- structure(list(ID = c(1L, 2L, 3L, 4L, 6L, 5L, 7L), dad = c(0L,
0L, 1L, 1L, 1L, 3L, 5L), mum = c(0L, 0L, 0L, 2L, 4L, 4L, 6L),
GEN = c(1L, 1L, 2L, 2L, 3L, 3L, 4L)), class = "data.frame", row.names = c(NA,
-7L))
hom<-function(data) {
library(Matrix)
library(foreach)
n<-max(as.numeric(fam[,"ID"]))
t<-min(as.numeric(fam[,"ID"]))
A<-Matrix(0,nrow=n,ncol=n, sparse=TRUE)
while(t <=n) {
s<-max(fam[t,"dad"],fam[t,"mum"])
d<-min(fam[t,"dad"],fam[t,"mum"])
if (s>0 & d>0 )
{
if (fam[t,"GEN"]==999 & s!=d)
{ warning("both dad and mum should be the same, different for at least one individual")
NULL
}
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)+0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]
foreach(j = 1:(t-1), .verbose=TRUE, .combine='c', .packages=c("Matrix", "foreach")) %do%
{
A[t,j]<- 0.5*(A[j,fam[t,"dad"]]+A[j,fam[t,"mum"]])
A[j,t]<- A[t,j]
}
}
if (s>0 & d==0 )
{
if ( fam[t,"GEN"]==999)
{ warning("both dad and mum should be the same, one parent equal to zero for at least individual")
NULL }
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
foreach(j = 1:(t-1), .verbose=TRUE, .combine='c', .packages=c("Matrix", "foreach")) %do%
{
A[t,j]<-0.5*A[j,s]
A[j,t]<-A[t,j]
}
}
if (s==0 )
{
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
}
cat(" MatbyGEN: ", t ,"\n")
t <- t+1
}
A
}
Output of the above example
%%MatrixMarket matrix coordinate real symmetric
7 7 26
1 1 1
3 1 .5
4 1 .5
5 1 .75
6 1 .5
7 1 .625
2 2 1
4 2 .5
5 2 .25
6 2 .25
7 2 .25
3 3 1.5
4 3 .25
5 3 .375
6 3 .875
7 3 .625
4 4 1.5
5 4 1
6 4 .875
7 4 .9375
5 5 1.8125
6 5 .6875
7 5 1.25
6 6 1.78125
7 6 1.234375
7 7 1.91796875
问题是让它更快地处理 300k x 300k 的矩阵，这需要几天或几周才能 运行 因为我已经 运行 宁它一段时间了，什么我可以做些什么让它 运行 更快吗？
N.B：将示例保存为"anything.txt"，然后将文件读取为"fam <- read.delim(, header = TRUE, sep="")"
你遇到的问题是这是递归的。每个循环都取决于前一个循环的结果。所以，你不能真正用向量化来解决问题。
如果你想为此使用 R，最好的办法是研究 Rcpp。我不太擅长 Rcpp，但我确实有一些建议。
最简单的做法是摆脱 foreach 循环并用常规 for 循环替换它。使用并行线程会产生很多开销，而且当函数是递归的时，工人很难靠自己真正做得更好。
# Before
foreach(j = 1:(t-1), .combine='c', .packages=c("Matrix", "foreach")) %do%
{ ... }
# After
for (j in 1:(t-1)) {
...
}
接下来要做的是考虑你是否真的需要一个稀疏矩阵。如果你没有内存问题，你还不如使用正则矩阵。
A<-Matrix(0,nrow=n,ncol=n, sparse=TRUE)
# to
A<-matrix(0,nrow=n,ncol=n)
最后要做的是重新考虑如何初始化所有内容。该代码的某些部分会重复多次，例如对 diag 的赋值。由于我们要对单独的元素求和，因此我们可以使用所有 3 个代码片段 2 - 0.5^(fam[t, 'GEN'] - 1) 共有的部分来初始化 diag。
A<-matrix(0,nrow=n,ncol=n)
diag(A) <- 2-0.5^(fam[["GEN"]]-1)
这很重要，因为这让我们可以跳过。您的原始代码片段有 1,000 行，其中 'mum' 和 'dad' 为 0。通过此初始化，我们可以直接跳到第一行，结果为 'mum' 或 'dad':
t_start <- min(which.max(fam$dad > 0), which.max(fam$mum > 0))
t_end <- max(fam[['ID']])
for (t in t_start:t_end) {
...
}
我决定跳过 if 语句，我想用 sum(c(..., ...)) 来总结一切。这样，如果子集的结果是 NULL，我仍然可以求和。总计：
t_start <- min(which.max(fam$dad > 0), which.max(fam$mum > 0))
t_end <- max(fam[['ID']])
A<-matrix(0,nrow=t_end,ncol=t_end)
diag(A) <- 2-0.5^(fam[["GEN"]]-1)
for (t in t_start:t_end) {
A[t,t]<- sum(c(A[t,t], 0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]))
for(j in 1:(t-1)) {
A[t,j]<- 0.5 * sum(c(A[j,fam[t,"dad"]],A[j,fam[t,"mum"]]))
A[j,t]<- A[t,j]
}
}
A
性能
Unit: microseconds
expr min lq mean median uq max neval
original 85759.901 86650.7515 88776.695 88740.050 90529.750 97433.2 100
non_foreach 47912.601 48528.5010 50699.867 50220.901 51782.651 88355.1 100
non_sparse_non_each 1423.701 1454.3015 1531.833 1471.451 1496.401 4126.3 100
final_change 953.102 981.8015 1212.264 1010.500 1026.052 21350.1 100
所有代码
fam <- structure(list(ID = c(1L, 2L, 3L, 4L, 6L, 5L, 7L), dad = c(0L,
0L, 1L, 1L, 1L, 3L, 5L), mum = c(0L, 0L, 0L, 2L, 4L, 4L, 6L),
GEN = c(1L, 1L, 2L, 2L, 3L, 3L, 4L)), class = "data.frame", row.names = c(NA,
-7L))
A<-matrix(0,nrow=7,ncol=7)
diag(A) <- 2-0.5^(fam[["GEN"]]-1)
t_start <- min(which.max(fam$dad > 0), which.max(fam$mum > 0))
t_end <- max(fam[['ID']])
for (t in t_start:t_end) {
A[t,t]<- sum(c(A[t,t], 0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]))
for(j in 1:(t-1)) {
A[t,j]<- 0.5 * sum(c(A[j,fam[t,"dad"]],A[j,fam[t,"mum"]]))
A[j,t]<- A[t,j]
}
}
A
hom<-function(data) {
library(Matrix)
library(foreach)
n<-max(as.numeric(fam[,"ID"]))
t<-min(as.numeric(fam[,"ID"]))
A<-Matrix(0,nrow=n,ncol=n, sparse=TRUE)
while(t <=n) {
s<-max(fam[t,"dad"],fam[t,"mum"])
d<-min(fam[t,"dad"],fam[t,"mum"])
if (s>0 & d>0 )
{
if (fam[t,"GEN"]==999 & s!=d)
{ warning("both dad and mum should be the same, different for at least one individual")
NULL
}
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)+0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]
foreach(j = 1:(t-1), .combine='c', .packages=c("Matrix", "foreach")) %do%
{
A[t,j]<- 0.5*(A[j,fam[t,"dad"]]+A[j,fam[t,"mum"]])
A[j,t]<- A[t,j]
}
}
if (s>0 & d==0 )
{
if ( fam[t,"GEN"]==999)
{ warning("both dad and mum should be the same, one parent equal to zero for at least individual")
NULL }
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
foreach(j = 1:(t-1), .combine='c', .packages=c("Matrix", "foreach")) %do%
{
A[t,j]<-0.5*A[j,s]
A[j,t]<-A[t,j]
}
}
if (s==0 )
{
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
}
# cat(" MatbyGEN: ", t ,"\n")
t <- t+1
}
A
}
hom2<-function(data) {
library(Matrix)
n<-max(as.numeric(fam[,"ID"]))
t<-min(as.numeric(fam[,"ID"]))
A<-Matrix(0,nrow=n,ncol=n, sparse = T)
while(t <=n) {
s<-max(fam[t,"dad"],fam[t,"mum"])
d<-min(fam[t,"dad"],fam[t,"mum"])
if (s>0 & d>0 )
{
if (fam[t,"GEN"]==999 & s!=d)
{ warning("both dad and mum should be the same, different for at least one individual")
NULL
}
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)+0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]
for (j in 1:(t-1)) {
A[t,j]<- 0.5*(A[j,fam[t,"dad"]]+A[j,fam[t,"mum"]])
A[j,t]<- A[t,j]
}
}
if (s>0 & d==0 )
{
if ( fam[t,"GEN"]==999)
{ warning("both dad and mum should be the same, one parent equal to zero for at least individual")
NULL }
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
for (j in 1:(t-1)) {
A[t,j]<-0.5*A[j,s]
A[j,t]<-A[t,j]
}
}
if (s==0 )
{
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
}
# cat(" MatbyGEN: ", t ,"\n")
t <- t+1
}
A
}
hom3<-function(data) {
n<-max(as.numeric(fam[,"ID"]))
t<-min(as.numeric(fam[,"ID"]))
A<-matrix(0,nrow=n,ncol=n)
while(t <=n) {
s<-max(fam[t,"dad"],fam[t,"mum"])
d<-min(fam[t,"dad"],fam[t,"mum"])
if (s>0 & d>0 )
{
if (fam[t,"GEN"]==999 & s!=d)
{ warning("both dad and mum should be the same, different for at least one individual")
NULL
}
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)+0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]
for (j in 1:(t-1)) {
A[t,j]<- 0.5*(A[j,fam[t,"dad"]]+A[j,fam[t,"mum"]])
A[j,t]<- A[t,j]
}
}
if (s>0 & d==0 )
{
if ( fam[t,"GEN"]==999)
{ warning("both dad and mum should be the same, one parent equal to zero for at least individual")
NULL }
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
for (j in 1:(t-1)) {
A[t,j]<-0.5*A[j,s]
A[j,t]<-A[t,j]
}
}
if (s==0 )
{
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
}
# cat(" MatbyGEN: ", t ,"\n")
t <- t+1
}
A
}
library(microbenchmark)
f_changed = function(fam) {
t_start <- min(which.max(fam$dad > 0), which.max(fam$mum > 0))
t_end <- max(fam[['ID']])
A<-matrix(0,nrow=t_end,ncol=t_end)
diag(A) <- 2-0.5^(fam[["GEN"]]-1)
for (t in t_start:t_end) {
A[t,t]<- sum(c(A[t,t], 0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]))
for(j in 1:(t-1)) {
A[t,j]<- 0.5 * sum(c(A[j,fam[t,"dad"]],A[j,fam[t,"mum"]]))
A[j,t]<- A[t,j]
}
}
A
}
microbenchmark(
original = {
hom(fam)
}
, non_foreach = {
hom2(fam)
}
, non_sparse_non_each = {
hom3(fam)
}
, final_change = {
f_changed(fam)
}
,times = 100
)
我正在尝试在 R 中计算一个 300,000x300,000 的矩阵，我的代码工作得很好，但现在已经 运行ning 好几天了，我怎样才能让它更高效、更省时？
我的代码运行良好，但已经 运行ning 好几天了，附件是我正在使用的代码的一个子集，ID 扩展到 300,000；我怎样才能使代码 运行 在几分钟内更快，因为它已经 运行 好几天了。
fam <- structure(list(ID = c(1L, 2L, 3L, 4L, 6L, 5L, 7L), dad = c(0L,
0L, 1L, 1L, 1L, 3L, 5L), mum = c(0L, 0L, 0L, 2L, 4L, 4L, 6L),
GEN = c(1L, 1L, 2L, 2L, 3L, 3L, 4L)), class = "data.frame", row.names = c(NA,
-7L))
hom<-function(data) {
library(Matrix)
library(foreach)
n<-max(as.numeric(fam[,"ID"]))
t<-min(as.numeric(fam[,"ID"]))
A<-Matrix(0,nrow=n,ncol=n, sparse=TRUE)
while(t <=n) {
s<-max(fam[t,"dad"],fam[t,"mum"])
d<-min(fam[t,"dad"],fam[t,"mum"])
if (s>0 & d>0 )
{
if (fam[t,"GEN"]==999 & s!=d)
{ warning("both dad and mum should be the same, different for at least one individual")
NULL
}
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)+0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]
foreach(j = 1:(t-1), .verbose=TRUE, .combine='c', .packages=c("Matrix", "foreach")) %do%
{
A[t,j]<- 0.5*(A[j,fam[t,"dad"]]+A[j,fam[t,"mum"]])
A[j,t]<- A[t,j]
}
}
if (s>0 & d==0 )
{
if ( fam[t,"GEN"]==999)
{ warning("both dad and mum should be the same, one parent equal to zero for at least individual")
NULL }
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
foreach(j = 1:(t-1), .verbose=TRUE, .combine='c', .packages=c("Matrix", "foreach")) %do%
{
A[t,j]<-0.5*A[j,s]
A[j,t]<-A[t,j]
}
}
if (s==0 )
{
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
}
cat(" MatbyGEN: ", t ,"\n")
t <- t+1
}
A
}
Output of the above example
%%MatrixMarket matrix coordinate real symmetric
7 7 26
1 1 1
3 1 .5
4 1 .5
5 1 .75
6 1 .5
7 1 .625
2 2 1
4 2 .5
5 2 .25
6 2 .25
7 2 .25
3 3 1.5
4 3 .25
5 3 .375
6 3 .875
7 3 .625
4 4 1.5
5 4 1
6 4 .875
7 4 .9375
5 5 1.8125
6 5 .6875
7 5 1.25
6 6 1.78125
7 6 1.234375
7 7 1.91796875
问题是让它更快地处理 300k x 300k 的矩阵，这需要几天或几周才能 运行 因为我已经 运行 宁它一段时间了，什么我可以做些什么让它 运行 更快吗？
N.B：将示例保存为"anything.txt"，然后将文件读取为"fam <- read.delim(, header = TRUE, sep="")"
你遇到的问题是这是递归的。每个循环都取决于前一个循环的结果。所以，你不能真正用向量化来解决问题。
如果你想为此使用 R，最好的办法是研究 Rcpp。我不太擅长 Rcpp，但我确实有一些建议。
最简单的做法是摆脱 foreach 循环并用常规 for 循环替换它。使用并行线程会产生很多开销，而且当函数是递归的时，工人很难靠自己真正做得更好。
# Before
foreach(j = 1:(t-1), .combine='c', .packages=c("Matrix", "foreach")) %do%
{ ... }
# After
for (j in 1:(t-1)) {
...
}
接下来要做的是考虑你是否真的需要一个稀疏矩阵。如果你没有内存问题，你还不如使用正则矩阵。
A<-Matrix(0,nrow=n,ncol=n, sparse=TRUE)
# to
A<-matrix(0,nrow=n,ncol=n)
最后要做的是重新考虑如何初始化所有内容。该代码的某些部分会重复多次，例如对 diag 的赋值。由于我们要对单独的元素求和，因此我们可以使用所有 3 个代码片段 2 - 0.5^(fam[t, 'GEN'] - 1) 共有的部分来初始化 diag。
A<-matrix(0,nrow=n,ncol=n)
diag(A) <- 2-0.5^(fam[["GEN"]]-1)
这很重要，因为这让我们可以跳过。您的原始代码片段有 1,000 行，其中 'mum' 和 'dad' 为 0。通过此初始化，我们可以直接跳到第一行，结果为 'mum' 或 'dad':
t_start <- min(which.max(fam$dad > 0), which.max(fam$mum > 0))
t_end <- max(fam[['ID']])
for (t in t_start:t_end) {
...
}
我决定跳过 if 语句，我想用 sum(c(..., ...)) 来总结一切。这样，如果子集的结果是 NULL，我仍然可以求和。总计：
t_start <- min(which.max(fam$dad > 0), which.max(fam$mum > 0))
t_end <- max(fam[['ID']])
A<-matrix(0,nrow=t_end,ncol=t_end)
diag(A) <- 2-0.5^(fam[["GEN"]]-1)
for (t in t_start:t_end) {
A[t,t]<- sum(c(A[t,t], 0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]))
for(j in 1:(t-1)) {
A[t,j]<- 0.5 * sum(c(A[j,fam[t,"dad"]],A[j,fam[t,"mum"]]))
A[j,t]<- A[t,j]
}
}
A
性能
Unit: microseconds
expr min lq mean median uq max neval
original 85759.901 86650.7515 88776.695 88740.050 90529.750 97433.2 100
non_foreach 47912.601 48528.5010 50699.867 50220.901 51782.651 88355.1 100
non_sparse_non_each 1423.701 1454.3015 1531.833 1471.451 1496.401 4126.3 100
final_change 953.102 981.8015 1212.264 1010.500 1026.052 21350.1 100
所有代码
fam <- structure(list(ID = c(1L, 2L, 3L, 4L, 6L, 5L, 7L), dad = c(0L,
0L, 1L, 1L, 1L, 3L, 5L), mum = c(0L, 0L, 0L, 2L, 4L, 4L, 6L),
GEN = c(1L, 1L, 2L, 2L, 3L, 3L, 4L)), class = "data.frame", row.names = c(NA,
-7L))
A<-matrix(0,nrow=7,ncol=7)
diag(A) <- 2-0.5^(fam[["GEN"]]-1)
t_start <- min(which.max(fam$dad > 0), which.max(fam$mum > 0))
t_end <- max(fam[['ID']])
for (t in t_start:t_end) {
A[t,t]<- sum(c(A[t,t], 0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]))
for(j in 1:(t-1)) {
A[t,j]<- 0.5 * sum(c(A[j,fam[t,"dad"]],A[j,fam[t,"mum"]]))
A[j,t]<- A[t,j]
}
}
A
hom<-function(data) {
library(Matrix)
library(foreach)
n<-max(as.numeric(fam[,"ID"]))
t<-min(as.numeric(fam[,"ID"]))
A<-Matrix(0,nrow=n,ncol=n, sparse=TRUE)
while(t <=n) {
s<-max(fam[t,"dad"],fam[t,"mum"])
d<-min(fam[t,"dad"],fam[t,"mum"])
if (s>0 & d>0 )
{
if (fam[t,"GEN"]==999 & s!=d)
{ warning("both dad and mum should be the same, different for at least one individual")
NULL
}
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)+0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]
foreach(j = 1:(t-1), .combine='c', .packages=c("Matrix", "foreach")) %do%
{
A[t,j]<- 0.5*(A[j,fam[t,"dad"]]+A[j,fam[t,"mum"]])
A[j,t]<- A[t,j]
}
}
if (s>0 & d==0 )
{
if ( fam[t,"GEN"]==999)
{ warning("both dad and mum should be the same, one parent equal to zero for at least individual")
NULL }
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
foreach(j = 1:(t-1), .combine='c', .packages=c("Matrix", "foreach")) %do%
{
A[t,j]<-0.5*A[j,s]
A[j,t]<-A[t,j]
}
}
if (s==0 )
{
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
}
# cat(" MatbyGEN: ", t ,"\n")
t <- t+1
}
A
}
hom2<-function(data) {
library(Matrix)
n<-max(as.numeric(fam[,"ID"]))
t<-min(as.numeric(fam[,"ID"]))
A<-Matrix(0,nrow=n,ncol=n, sparse = T)
while(t <=n) {
s<-max(fam[t,"dad"],fam[t,"mum"])
d<-min(fam[t,"dad"],fam[t,"mum"])
if (s>0 & d>0 )
{
if (fam[t,"GEN"]==999 & s!=d)
{ warning("both dad and mum should be the same, different for at least one individual")
NULL
}
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)+0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]
for (j in 1:(t-1)) {
A[t,j]<- 0.5*(A[j,fam[t,"dad"]]+A[j,fam[t,"mum"]])
A[j,t]<- A[t,j]
}
}
if (s>0 & d==0 )
{
if ( fam[t,"GEN"]==999)
{ warning("both dad and mum should be the same, one parent equal to zero for at least individual")
NULL }
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
for (j in 1:(t-1)) {
A[t,j]<-0.5*A[j,s]
A[j,t]<-A[t,j]
}
}
if (s==0 )
{
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
}
# cat(" MatbyGEN: ", t ,"\n")
t <- t+1
}
A
}
hom3<-function(data) {
n<-max(as.numeric(fam[,"ID"]))
t<-min(as.numeric(fam[,"ID"]))
A<-matrix(0,nrow=n,ncol=n)
while(t <=n) {
s<-max(fam[t,"dad"],fam[t,"mum"])
d<-min(fam[t,"dad"],fam[t,"mum"])
if (s>0 & d>0 )
{
if (fam[t,"GEN"]==999 & s!=d)
{ warning("both dad and mum should be the same, different for at least one individual")
NULL
}
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)+0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]
for (j in 1:(t-1)) {
A[t,j]<- 0.5*(A[j,fam[t,"dad"]]+A[j,fam[t,"mum"]])
A[j,t]<- A[t,j]
}
}
if (s>0 & d==0 )
{
if ( fam[t,"GEN"]==999)
{ warning("both dad and mum should be the same, one parent equal to zero for at least individual")
NULL }
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
for (j in 1:(t-1)) {
A[t,j]<-0.5*A[j,s]
A[j,t]<-A[t,j]
}
}
if (s==0 )
{
A[t,t]<- 2-0.5^(fam[t,"GEN"]-1)
}
# cat(" MatbyGEN: ", t ,"\n")
t <- t+1
}
A
}
library(microbenchmark)
f_changed = function(fam) {
t_start <- min(which.max(fam$dad > 0), which.max(fam$mum > 0))
t_end <- max(fam[['ID']])
A<-matrix(0,nrow=t_end,ncol=t_end)
diag(A) <- 2-0.5^(fam[["GEN"]]-1)
for (t in t_start:t_end) {
A[t,t]<- sum(c(A[t,t], 0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]))
for(j in 1:(t-1)) {
A[t,j]<- 0.5 * sum(c(A[j,fam[t,"dad"]],A[j,fam[t,"mum"]]))
A[j,t]<- A[t,j]
}
}
A
}
microbenchmark(
original = {
hom(fam)
}
, non_foreach = {
hom2(fam)
}
, non_sparse_non_each = {
hom3(fam)
}
, final_change = {
f_changed(fam)
}
,times = 100
)