From 7d6718046c2a8f3276392bfd35fcafc02d0296a7 Mon Sep 17 00:00:00 2001
From: Pierre Santagostini <pierre.santagostini@agrocampus-ouest.fr>
Date: Thu, 4 Jul 2024 20:08:16 +0200
Subject: [PATCH] Use of data.table::rbindlist() instead of rbind()
---
R/lauricella.R | 155 +++++++++++++++++++++++++++++++++----------------
1 file changed, 104 insertions(+), 51 deletions(-)
diff --git a/R/lauricella.R b/R/lauricella.R
index 541329e..bb5dca0 100644
--- a/R/lauricella.R
+++ b/R/lauricella.R
@@ -34,8 +34,9 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
#' IEEE Signal Processing Letters, vol. 26 no. 7, July 2019.
#' \doi{10.1109/LSP.2019.2915000}
#' @importFrom utils combn
+ #' @importFrom data.table rbindlist
#' @export
-
+
# Number of variables
n <- length(x)
@@ -56,7 +57,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# M: data.frame of the indices for the nested sums
# (i.e. all arrangements of n elements from {0:k})
M <- expand.grid(rep(list(0:k), n))
-
+
# Sum of the indices
Msum <- rowSums(M)
@@ -100,7 +101,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
lnapoch <- sapply(Msumunique, function(j) lnpochhammer(a, j))
names(lnapoch) <- Msumunique
-
+
# Logarithms of pochhammer(b_i, m_i) for m_1 = 0...k, ..., m_n = 0...k
lnbpoch <- as.data.frame(
matrix(nrow = max(M) + 1, ncol = n, dimnames = list(Munique, 1:n)))
@@ -163,39 +164,52 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
k1 <- k + (1:kstep)
result <- result + res2
- # # M: data.frame of the indices for the nested sums
+ buildMlist <- function(j, isupp, k = k, k1 = k1) {
+ jsupp <- isupp[, j]
+ Ml <- rep(list(0:k), n)
+ for (j in jsupp)
+ Ml[[j]] <- k1
+ return(expand.grid(Ml))
+ }
+
+ buildxlist <- function(j, isupp, x, xsupp) {
+ jsupp <- isupp[, j]
+ xl <- rep(list(x), n)
+ xl[jsupp] <- xsupp[jsupp]
+ xl[-jsupp] <- x[-jsupp]
+ return(expand.grid(xl))
+ }
+
+ # M: data.frame of the indices for the nested sums
# M <- expand.grid(rep(list(k1), n))
# if (n > 1) {
# for (i in 1:(n-1)) {
- # Mlist <- c( rep(list(0:k), n-i), rep(list(k1), i) )
- # M <- rbind( M, expand.grid(Mlist) )
- # for (j in 1:(n-1)) {
- # Mlist <- Mlist[c(n, 1:(n-1))]
+ # indsupp <- combn(n, i)
+ # for (j in 1:ncol(indsupp)) {
+ # jsupp <- indsupp[, j]
+ # Mlist <- vector("list", n)
+ # for (l in jsupp) Mlist[[l]] <- k1
+ # for (l in (1:n)[-jsupp]) Mlist[[l]] <- 0:k
# M <- rbind(M, expand.grid(Mlist))
# }
# }
# }
- # M1 <- M
-
- # M: data.frame of the indices for the nested sums
- M <- expand.grid(rep(list(k1), n))
+ Mlist <- list(expand.grid(rep(list(k1), n)))
+ if (n == 1) {
+ M <- Mlist[[1]]
+ }
if (n > 1) {
for (i in 1:(n-1)) {
indsupp <- combn(n, i)
- for (j in 1:ncol(indsupp)) {
- jsupp <- indsupp[, j]
- Mlist <- vector("list", n)
- for (l in jsupp) Mlist[[l]] <- k1
- for (l in (1:n)[-jsupp]) Mlist[[l]] <- 0:k
- M <- rbind(M, expand.grid(Mlist))
- }
+ Mlist <- c(Mlist,
+ lapply(1:ncol(indsupp), buildMlist, isupp = indsupp, k = k, k1 = k1))
}
+ M <- data.frame(rbindlist(Mlist))
}
- # M2 <- M
# Sum of the indices
Msum <- rowSums(M)
-
+
Munique <- (max(Munique)+1):max(M)
Msumunique <- (max(Msumunique)+1):max(Msum)
@@ -208,17 +222,30 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
xfactsupp[i, j] <- x[j]^Munique[i]
}
- gridxfact <- expand.grid(xfactsupp)
- for (i in 1:(n-1)) {
- indsupp <- combn(n, i)
- for (j in 1:ncol(indsupp)) {
- jsupp <- indsupp[, j]
- xfactlist <- vector("list", n)
- names(xfactlist) <- names(gridxfact)
- xfactlist[jsupp] <- xfactsupp[jsupp]
- xfactlist[-jsupp] <- xfact[-jsupp]
- gridxfact <- rbind(gridxfact, expand.grid(xfactlist))
+ # gridxfact <- expand.grid(xfactsupp)
+ # for (i in 1:(n-1)) {
+ # indsupp <- combn(n, i)
+ # for (j in 1:ncol(indsupp)) {
+ # jsupp <- indsupp[, j]
+ # xfactlist <- vector("list", n)
+ # names(xfactlist) <- names(gridxfact)
+ # xfactlist[jsupp] <- xfactsupp[jsupp]
+ # xfactlist[-jsupp] <- xfact[-jsupp]
+ # gridxfact <- rbind(gridxfact, expand.grid(xfactlist))
+ # }
+ # }
+ xfactlist <- list(expand.grid(xfactsupp))
+ if (n == 1) {
+ gridxfact <- xfactlist[[1]]
+ }
+ if (n > 1) {
+ for (i in 1:(n-1)) {
+ indsupp <- combn(n, i)
+ xfactlist <- c(xfactlist,
+ lapply(1:ncol(indsupp), buildxlist, isupp = indsupp, x = xfact, xsupp = xfactsupp))
}
+ names(xfactlist[[1]]) <- names(xfactlist[[2]])
+ gridxfact <- data.frame(rbindlist(xfactlist))
}
prodxfact <- apply(gridxfact, 1, prod)
@@ -227,17 +254,30 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
lnfactsupp <- as.data.frame(
matrix(lfactorial(Munique), nrow = length(Munique),
ncol = n, dimnames = list(Munique, 1:n)))
- gridlnfact <- expand.grid(lnfactsupp)
- for (i in 1:(n-1)) {
- indsupp <- combn(n, i)
- for (j in 1:ncol(indsupp)) {
- jsupp <- indsupp[, j]
- lnfactlist <- vector("list", n)
- names(lnfactlist) <- names(gridlnfact)
- lnfactlist[jsupp] <- lnfactsupp[jsupp]
- lnfactlist[-jsupp] <- lnfact[-jsupp]
- gridlnfact <- rbind(gridlnfact, expand.grid(lnfactlist))
+ # gridlnfact <- expand.grid(lnfactsupp)
+ # for (i in 1:(n-1)) {
+ # indsupp <- combn(n, i)
+ # for (j in 1:ncol(indsupp)) {
+ # jsupp <- indsupp[, j]
+ # lnfactlist <- vector("list", n)
+ # names(lnfactlist) <- names(gridlnfact)
+ # lnfactlist[jsupp] <- lnfactsupp[jsupp]
+ # lnfactlist[-jsupp] <- lnfact[-jsupp]
+ # gridlnfact <- rbind(gridlnfact, expand.grid(lnfactlist))
+ # }
+ # }
+ lnfactlist <- list(expand.grid(lnfactsupp))
+ if (n == 1) {
+ gridlnfact <- lnfactlist[[1]]
+ }
+ if (n > 1) {
+ for (i in 1:(n-1)) {
+ indsupp <- combn(n, i)
+ lnfactlist <- c(lnfactlist,
+ lapply(1:ncol(indsupp), buildxlist, isupp = indsupp, x = lnfact, xsupp = lnfactsupp))
}
+ names(lnfactlist[[1]]) <- names(lnfactlist[[2]])
+ gridlnfact <- data.frame(rbindlist(lnfactlist))
}
sumlnfact <- rowSums(gridlnfact)
@@ -260,17 +300,30 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
lnbpochsupp[i, j] <- lnpochhammer(b[j], Munique[i])
}
# Table of all combinations of the log(pochhammer(b_i, m_i))
- gridlnbpoch <- expand.grid(lnbpochsupp)
- for (i in 1:(n-1)) {
- indsupp <- combn(n, i)
- for (j in 1:ncol(indsupp)) {
- jsupp <- indsupp[, j]
- lnbpochlist <- vector("list", n)
- names(lnbpochlist) <- names(gridlnbpoch)
- lnbpochlist[jsupp] <- lnbpochsupp[jsupp]
- lnbpochlist[-jsupp] <- lnbpoch[-jsupp]
- gridlnbpoch <- rbind(gridlnbpoch, expand.grid(lnbpochlist))
+ # gridlnbpoch <- expand.grid(lnbpochsupp)
+ # for (i in 1:(n-1)) {
+ # indsupp <- combn(n, i)
+ # for (j in 1:ncol(indsupp)) {
+ # jsupp <- indsupp[, j]
+ # lnbpochlist <- vector("list", n)
+ # names(lnbpochlist) <- names(gridlnbpoch)
+ # lnbpochlist[jsupp] <- lnbpochsupp[jsupp]
+ # lnbpochlist[-jsupp] <- lnbpoch[-jsupp]
+ # gridlnbpoch <- rbind(gridlnbpoch, expand.grid(lnbpochlist))
+ # }
+ # }
+ lnbpochlist <- list(expand.grid(lnbpochsupp))
+ if (n == 1) {
+ gridlnbpoch <- lnbpochlist[[1]]
+ }
+ if (n > 1) {
+ for (i in 1:(n-1)) {
+ indsupp <- combn(n, i)
+ lnbpochlist <- c(lnbpochlist,
+ lapply(1:ncol(indsupp), buildxlist, isupp = indsupp, x = lnbpoch, xsupp = lnbpochsupp))
}
+ names(lnbpochlist[[1]]) <- names(lnbpochlist[[2]])
+ gridlnbpoch <- data.frame(rbindlist(lnbpochlist))
}
# Sum of the logarithms
sumlnpochb <- rowSums(gridlnbpoch)
--
GitLab