
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.gh <- volta.quad(md.cov.ck,n.pontos=7,Variancia=FALSE) preditos.gh <- data.frame(preditos.gh) write.table(preditos.gh,"pred_by4k7ns1000MBM.txt") 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. ------------------------------------------------------------------- mat.cokri <- md.cov.ck 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. Meu email é beleti.junior@gmail.com 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.