Olá Benilton,
Primeiramente gostaria de agradecer novamente pela paciência e pela atenção, e me desculpar por não ter conhecimento necessário sobre funções e procedimentos estatísticos.
O foco da minha pesquisa é o processamento paralelo. Estou utilizando scripts, (tais como de obtenção dos valores preditos utilizando os dados do pivo), pois foi constatado problemas quanto ao tempo de processamento de algumas rotinas desses scripts.
Na rotina de obtenção dos valores preditos utilizando os dados do pivo (trecho a seguir), o problema quanto a tempo de processamento está na chamada "preditos.simu <- volta.cokri(
md.cov.ck,num.simu=1000,int.conf=0.95)". (lembrando que tais scripts necessitam do pacote geoComp)
-------------------------------------------------------------------
require(geoComp)
require(MASS)
require(statmod)
require(geoR)
data(pivo)
dados <- pivo[,c(6,7,8,1,2)]
dados <- as.geoComp(dados)
bor <- cbind(c(0,seq(0,200,l=100),0),c(0,sqrt(200^2-seq(0,200,l=100)^2),0))
estima <- mec(dados)
gr <- pred_grid(bor, by=4)
md.cov.ck <- cokrigagem(estima[[1]]$Estimativas, loc=gr,dados.comp=dados)
preditos.simu <- volta.cokri(
md.cov.ck,num.simu=1000,int.conf=0.95) #Ponto com problema quanto a processamento
preditos.simu <- data.frame(preditos.simu[[1]])
preditos.simu.ic <- data.frame(preditos.simu[[2]])
write.table(preditos.simu,"predsimu_by4k7ns1000MBM.txt")
write.table(preditos.simu.ic,"predsimuic_by4k7ns1000MBM.txt")
-------------------------------------------------------------------
Visando localizar o ponto em que a execução gasta um maior tempo de processamento, "entrei" na função volta.cokri(código a seguir) e fui investigar a causa da "demora" na execução dessa função.
-------------------------------------------------------------------
volta.cokri <-
function(mat.cokri, num.simu, retorna.tudo=FALSE, int.conf=0.95){
nlinhas <- dim(mat.cokri[[1]])[1]
compos1 <- data.frame(matrix(nrow=nlinhas/2,ncol=3))
compos <- data.frame(matrix(nrow=nlinhas/2,ncol=3))
for(i in 1:num.simu){
g <- mvrnorm(n=1, mat.cokri[[1]],mat.cokri[[2]])
seq1 <- seq(1,nlinhas,by=2)
seq2 <- seq(2,nlinhas,by=2)
y1 <- g[seq1]
y2 <- g[seq2]
gerado <- data.frame(y1,y2)
compos <- agl(gerado)
compos1 <- cbind(compos,compos1)
}
compos2 <- as.matrix(compos1)
dim.vetor <- num.simu*3
sy1 <- seq(1,dim.vetor,by=3)
sy2 <- seq(2,dim.vetor,by=3)
sy3 <- seq(3,dim.vetor,by=3)
amostra1 <- as.matrix(compos2[,sy1],ncol=num.simu)
amostra2 <- as.matrix(compos2[,sy2],ncol=num.simu)
amostra3 <- as.matrix(compos2[,sy3],ncol=num.simu)
if(retorna.tudo == TRUE){
retorna <- list()
retorna[[1]] <- amostra1
retorna[[2]] <- amostra2
retorna[[3]] <- amostra3
return(retorna)
}
if(retorna.tudo == FALSE){
med1 <- apply(amostra1,1,mean)
med2 <- apply(amostra2,1,mean)
med3 <- apply(amostra3,1,mean)
q1 <- t(apply(amostra1,1,quantile,prob=c(1-int.conf,int.conf)))
q2 <- t(apply(amostra2,1,quantile,prob=c(1-int.conf,int.conf)))
q3 <- t(apply(amostra3,1,quantile,prob=c(1-int.conf,int.conf)))
quantis <- cbind(q1,q2,q3)
resultado <- list()
resultado$preditos <- data.frame(med1,med2,med3)
names(resultado$preditos) <- c("Areia","Silte","Argila")
resultado$intervalo <- data.frame(quantis)
names(resultado$intervalo) <- c("LI Areia", "LS Areia", "LI Silte", "LS Silte", "LI Argila", "LS Argila")
return(resultado)}
}
-------------------------------------------------------------------
Na função volta.cokri verifiquei que o elevado tempo de processamento se concentrava na chamada "g <- mvrnorm(n=1, mat.cokri[[1]],mat.cokri[[2]])". Então desmembrei a função para realizar experimentos com pacotes paralelos afim de criar modelos de execução paralelos.
A seguir está o trecho da função desmembrada.
-------------------------------------------------------------------
retorna.tudo=FALSE
int.conf=0.95
nlinhas <- dim(mat.cokri[[1]])[1]
compos1 <- data.frame(matrix(nrow=nlinhas/2,ncol=3))
compos <- data.frame(matrix(nrow=nlinhas/2,ncol=3))
mvrnorm <- function(n = 1, mu, Sigma, tol=1e-6, empirical = FALSE)
{
p <- length(mu)
if(!all(dim(Sigma) == c(p,p))) stop("incompatible arguments")
eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
ev <- eS$values
if(!all(ev >= -tol*abs(ev[1L]))) ev <- ev *(-1)
X <- matrix(rnorm(p * n), n)
if(empirical) {
X <- scale(X, TRUE, FALSE)
X <- X %*% svd(X, nu = 0)$v
X <- scale(X, FALSE, TRUE)
}
X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X)
nm <- names(mu)
if(is.null(nm) && !is.null(dn <- dimnames(Sigma))) nm <- dn[[1L]]
dimnames(X) <- list(nm, NULL)
if(n == 1) drop(X) else t(X)
}
m1 <- mat.cokri[[2]][1:2601,1:2601] #divisao da matriz original
m2 <- mat.cokri[[2]][2602:5202,1:2601]
m3 <- mat.cokri[[2]][1:2601,2602:5202]
m4 <- mat.cokri[[2]][2602:5202,2602:5202]
#chamada a funcao mclapply (neste exemplo esta com 1:4, mas que na execução real seria 1:1000)
mclapply(1:4, function(i) {
g1 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m1)
g2 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m2)
g3 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m3)
g4 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m4)
list(g1, g2, g3, g4)
},mc.cores=4
)
g <- list(list(1:2601, g1), list(1:2601, g2), list(1:2601, g3),list(1:2601, g4))
seq1 <- seq(1,nlinhas,by=2)
seq2 <- seq(2,nlinhas,by=2)
y1 <- g[seq1]
y2 <- g[seq2]
gerado <- data.frame(y1,y2)
compos <- agl(gerado)
compos1 <- cbind(compos,compos1)
compos2 <- as.matrix(compos1)
-------------------------------------------------------------------
Esses trechos de código podem reproduzir o problema. Se não tiver o pacote geoComp, posso lhe enviar por email.
Sobre a possivel solução apresentada, consegui aplica-lá.
Mais uma vez agradeço a sua disposição e peço desculpas pelo encomodo.
Att, Carlos Beleti.