From bd49fe4dbfbf316767fa5dd169e4ca42636c54b7 Mon Sep 17 00:00:00 2001
From: Pierre Santagostini <pierre.santagostini@agrocampus-ouest.fr>
Date: Mon, 8 Jul 2024 19:25:20 +0200
Subject: [PATCH] Modifications to reduce the execution time.
---
R/dlauricella.R | 181 +++++++++++++++++++++++++++---------------------
1 file changed, 102 insertions(+), 79 deletions(-)
diff --git a/R/dlauricella.R b/R/dlauricella.R
index ede99d1..dcb6784 100644
--- a/R/dlauricella.R
+++ b/R/dlauricella.R
@@ -1,9 +1,7 @@
-#' @importFrom data.table rbindlist
-
dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
-
+
p <- length(lambda)
-
+
buildMlist <- function(j, isupp, k, k1, p = p) {
jsupp <- isupp[, j]
Ml <- rep(list(0:k), p)
@@ -11,7 +9,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
Ml[[j]] <- k1
return(expand.grid(Ml))
}
-
+
buildxlist <- function(j, isupp, x, xsupp, p = p) {
jsupp <- isupp[, j]
xl <- rep(list(x), p)
@@ -19,29 +17,29 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
xl[-jsupp] <- x[-jsupp]
return(expand.grid(xl))
}
-
+
prodlambda <- prod(lambda)
-
+
lambdanu <- lambda*nu1/nu2
prodlambdanu <- prod(lambdanu)
-
+
k <- 5
-
+
# 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), p))
M <- M[-1, , drop = FALSE]
-
+
Munique <- 0:k
-
+
# Sum of the indices
Msum <- rowSums(M)
- # Msumunique <- 0:max(Msum)
-
+ Msumunique <- 0:max(Msum)
+
kstep <- 5
-
+
if (lambdanu[p] <= 1) { # lambda[1] < ... < lambda[p] < 1
-
+
if (p > 1 & lambda[p] == 1) {
lambda <- lambda[1:(p-1)]
p1 <- p-1
@@ -52,25 +50,25 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
} else {
p1 <- p
}
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- matrix(rep(Mpoch, p1), ncol = p1, byrow = FALSE)
matM <- matrix(rep(Munique, p1), ncol = p1, byrow = FALSE)
matlabmda <- matrix(rep(log(1 - lambdanu[1:p1]), length(Munique)), ncol = p1, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p1), ncol = p1, byrow = FALSE)
-
+
# matcommun <- exp(matpoch + lambdaM - matfact)
matcommun <- as.data.frame(matpoch + lambdaM - matfact)
# matcommun <- as.data.frame(apply(matcommun, 2, Re))
matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
gridcommun <- expand.grid(matcommun)
gridcommun <- gridcommun[-1, , drop = FALSE]
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- sum(
# commun * sapply(Msum, function(i) {
# pochhammer(1, i) / ( pochhammer((nu1 + p)/2, i) * i )
@@ -81,7 +79,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
})
)))
-
+
# Next elements of the sum, until the expected precision
k1 <- 1:k
derive <- 0
@@ -90,7 +88,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
k <- k1[length(k1)]
k1 <- k + (1:kstep)
derive <- derive + d
-
+
# M: data.frame of the indices for the nested sums
# M <- expand.grid(rep(list(k1), p1))
# if (p1 > 1) {
@@ -117,12 +115,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
M <- data.frame(rbindlist(Mlist))
}
-
+
# Sum of the indices
Msum <- rowSums(M)
-
+
Munique <- (max(Munique)+1):max(M)
- # Msumunique <- (max(Msumunique)+1):max(Msum)
+ Msumunique <- unique(Msum)
+ names(Msumunique) <- as.character(Msumunique)
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- matrix(rep(Mpoch, p1), ncol = p1, byrow = FALSE)
@@ -130,13 +129,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
matlabmda <- matrix(rep(log(1 - lambdanu[1:p1]), length(Munique)), ncol = p1, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p1), ncol = p1, byrow = FALSE)
-
+
matcommun <- rbind(matcommun, matcommunsupp)
# matcommunsupp <- exp(matpoch + lambdaM - matfact)
matcommunsupp <- as.data.frame(matpoch + lambdaM - matfact)
# matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
communlist <- list(expand.grid(matcommunsupp))
if (p1 == 1) {
gridcommun <- communlist[[1]]
@@ -152,49 +151,53 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- sum(
# commun * sapply(Msum, function(i) {
# pochhammer(1, i) / ( pochhammer((nu1 + p)/2, i) * i )
# })
# )
- d <- Re(sum(exp(
- commun + sapply(Msum, function(i) {
- lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
- })
- )))
-
+ # d <- Re(sum(exp(
+ # commun + sapply(Msum, function(i) {
+ # lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
+ # })
+ # )))
+ pochMsum <- sapply(Msumunique, function(i) {
+ lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
+ })
+ names(pochMsum) <- as.character(Msumunique)
+ d <- Re(sum(exp(commun + pochMsum[as.character(Msum)])))
}
-
+
derive <- as.numeric(-0.5 * log(prodlambda) - (nu2 + p)/2 * derive)
-
+
} else if (lambdanu[1] > 1) { # 1 < lambda[1] < ... < lambda[p]
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- matrix(rep(Mpoch, p), ncol = p, byrow = FALSE)
matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
matlabmda <- matrix(rep(log(1 - 1/lambdanu), length(Munique)), ncol = p, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
-
+
# matcommun <- exp(matpoch + lambdaM - matfact)
matcommun <- as.data.frame(matpoch + lambdaM - matfact)
# matcommun <- as.data.frame(apply(matcommun, 2, Re))
matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
gridcommun <- expand.grid(matcommun)
gridcommun <- gridcommun[-1, , drop = FALSE]
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
d <- 0
for (i in 1:length(Msum)) {
A <- sum(1/(0:(Msum[i]-1) + (nu1+p)/2))
d <- d - exp(commun[i]) * A
}
d <- Re(d)
-
+
# Next elements of the sum, until the expected precision
k1 <- 1:k
derive <- 0
@@ -205,7 +208,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
k1 <- k + (1:kstep)
derive <- derive + d
# vd <- c(vd, d); vderive <- c(vderive, derive)
-
+
# # M: data.frame of the indices for the nested sums
# M <- as.data.frame(matrix(nrow = 0, ncol = p))
# if (p > 1) {
@@ -219,7 +222,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# }
# }
# M <- rbind( M, expand.grid(rep(list(k1), p)) )
-
+
# M: data.frame of the indices for the nested sums
# M <- expand.grid(rep(list(k1), p))
# if (p > 1) {
@@ -246,12 +249,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
M <- data.frame(rbindlist(Mlist))
}
-
+
Msum <- rowSums(M)
-
+
Munique <- (max(Munique) + 1):max(M)
- # Msumunique <- (max(Msumunique) + 1):max(Msum)
-
+ Msumunique <- unique(Msum)
+ names(Msumunique) <- as.character(Msumunique)
+
# d <- 0
# for (i in 1:length(Msum)) {
# commun <- prod(
@@ -262,20 +266,20 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# A <- sum(1/(0:(Msum[i]-1) + (1+p)/2))
# d <- d - commun * A # / pochhammer((1 + p)/2, Msum[i])
# }
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- matrix(rep(Mpoch, p), ncol = p, byrow = FALSE)
matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
matlabmda <- matrix(rep(log(1 - 1/lambdanu), length(Munique)), ncol = p, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
-
+
matcommun <- rbind(matcommun, matcommunsupp)
# matcommunsupp <- exp(matpoch + lambdaM - matfact)
matcommunsupp <- as.data.frame(matpoch + lambdaM - matfact)
# matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
communlist <- list(expand.grid(matcommunsupp))
if (p == 1) {
gridcommun <- communlist[[1]]
@@ -291,22 +295,29 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- 0
# for (i in 1:length(Msum)) {
# A <- sum(1/(0:(Msum[i]-1) + (1+p)/2))
# d <- d - commun[i] * A
# }
- A <- sapply(Msum, function (i) sum(1/(0:(i-1) + (nu1+p)/2)))
+ # d <- commun + pochMsum[Msum]
+ # A <- sapply(Msum, function (i) sum(1/(0:(i-1) + (nu1+p)/2)))
+ pochMsum <- sapply(Msumunique, function(i) {
+ sum(1/(0:(i-1) + (nu1+p)/2))
+ })
+ names(pochMsum) <- as.character(Msumunique)
+ A <- pochMsum[as.character(Msum)]
+
d <- Re(-sum(exp(commun)*A))
-
+
}
-
+
derive <- as.numeric(-0.5 * log(prodlambda) - (nu2 + p)/2 * prod(1/sqrt(lambdanu)) * derive)
-
+
} else { # lambda[1] < ... < 1 < ... < lambda[p]
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- cbind(matrix(rep(Mpoch, p-1), ncol = p-1, byrow = FALSE), 0)
# Mpochp <- sapply(Munique, function(j) lnpochhammer(nu1/2, j))
@@ -315,18 +326,18 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
matlabmda <- matrix(rep(log(1 - c(lambdanu[-p], 1)/lambdanu[p]), length(Munique)), ncol = p, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
-
+
# matcommun <- exp(matpoch + lambdaM - matfact)
matcommun <- as.data.frame(matpoch + lambdaM - matfact)
# matcommun <- as.data.frame(apply(matcommun, 2, Re))
matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
gridcommun <- expand.grid(matcommun)
gridcommun <- gridcommun[-1, , drop = FALSE]
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- sum(
# commun * sapply(1:length(Msum), function(i) {
# pochhammer(0.5, M[i, p]) * pochhammer(1, Msum[i]) / ( pochhammer((nu1 + p)/2, Msum[i]) * Msum[i] )
@@ -338,7 +349,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
( lnpochhammer((nu1 + p)/2, Msum[i]) + log(Msum[i]) )
})
)))
-
+
# Next elements of the sum, until the expected precision
k1 <- 1:k
derive <- 0
@@ -347,7 +358,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
k <- k1[length(k1)]
k1 <- k + (1:kstep)
derive <- derive + d
-
+
# # M: data.frame of the indices for the nested sums
# M <- as.data.frame(matrix(nrow = 0, ncol = p))
# if (p > 1) {
@@ -361,7 +372,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# }
# }
# M <- rbind( M, expand.grid(rep(list(k1), p)) )
-
+
# M <- expand.grid(rep(list(k1), p))
# if (p > 1) {
# for (i in 1:(p-1)) {
@@ -375,7 +386,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# }
# }
# }
-
+
# M: data.frame of the indices for the nested sums
Mlist <- list(expand.grid(rep(list(k1), p)))
if (p == 1) {
@@ -389,12 +400,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
M <- data.frame(rbindlist(Mlist))
}
-
+
Msum <- rowSums(M)
-
+
Munique <- (max(Munique)+1):max(M)
- # Msumunique <- (max(Msumunique)+1):max(Msum)
-
+ Msumunique <- unique(Msum)
+ names(Msumunique) <- as.character(Msumunique)
+
# d <- 0
# for (i in 1:length(Msum)) {
# commun <- prod(
@@ -405,7 +417,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# commun <- commun*(1 - 1/lambda[p])^M[i, p]/factorial(M[i, p])
# d <- d + commun * pochhammer(0.5, M[i, p])*pochhammer(1, Msum[i]) / ( pochhammer((1 + p)/2, Msum[i]) * Msum[i] )
# }
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- cbind(matrix(rep(Mpoch, p-1), ncol = p-1, byrow = FALSE), 0)
# Mpochp <- sapply(Munique, function(j) lnpochhammer(nu1/2, j))
@@ -414,13 +426,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
matlabmda <- matrix(rep(log(1 - c(lambdanu[-p], 1)/lambdanu[p]), length(Munique)), ncol = p, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
-
+
matcommun <- rbind(matcommun, matcommunsupp)
# matcommunsupp <- exp(matpoch + lambdaM - matfact)
matcommunsupp <- as.data.frame(matpoch + lambdaM - matfact)
# matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
communlist <- list(expand.grid(matcommunsupp))
if (p == 1) {
gridcommun <- communlist[[1]]
@@ -436,27 +448,38 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- sum(
# commun * sapply(1:length(Msum), function(i) {
# pochhammer(0.5, M[i, p]) * pochhammer(1, Msum[i]) / ( pochhammer((nu1 + p)/2, Msum[i]) * Msum[i] )
# })
# )
+ # d <- Re(sum(exp(
+ # commun + sapply(1:length(Msum), function(i) {
+ # lnpochhammer(nu1/2, M[i, p]) + lnpochhammer(1, Msum[i]) -
+ # ( lnpochhammer((nu1 + p)/2, Msum[i]) + log(Msum[i]) )
+ # })
+ # )))
+ pochnup <- sapply(0:max(M), function(i) {
+ lnpochhammer(nu1/2, i)
+ })
+ names(pochnup) <- as.character(0:max(M))
+ pochMsum <- sapply(Msumunique, function(i) {
+ lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
+ })
+ names(pochMsum) <- as.character(Msumunique)
d <- Re(sum(exp(
- commun + sapply(1:length(Msum), function(i) {
- lnpochhammer(nu1/2, M[i, p]) + lnpochhammer(1, Msum[i]) -
- ( lnpochhammer((nu1 + p)/2, Msum[i]) + log(Msum[i]) )
- })
+ commun + pochMsum[as.character(Msum)] + pochnup[as.character(M[, p])]
)))
-
+
}
-
+
derive <- as.numeric(-0.5 * log(prodlambda) + (nu2 + p)/2 * (log(lambdanu[p]) - derive))
-
+
}
-
+
attr(derive, "epsilon") <- eps
attr(derive, "k") <- k
-
+
return(derive)
}
--
GitLab