<div>Olá Benilton,</div><div><br></div><div>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. </div>
<div><br></div><div>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.</div>
<div><br></div><div>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(<a href="http://md.cov.ck">md.cov.ck</a>,num.simu=1000,int.conf=0.95)". (lembrando que tais scripts necessitam do pacote geoComp)</div>
<div><br></div><div>-------------------------------------------------------------------</div><div>require(geoComp)</div><div>require(MASS)</div><div>require(statmod)</div><div>require(geoR)</div><div>data(pivo)</div><div>
dados <- pivo[,c(6,7,8,1,2)]</div><div>dados <- as.geoComp(dados)</div><div>bor <- cbind(c(0,seq(0,200,l=100),0),c(0,sqrt(200^2-seq(0,200,l=100)^2),0)) </div><div>estima <- mec(dados)</div><div>gr <- pred_grid(bor, by=4)</div>
<div><a href="http://md.cov.ck">md.cov.ck</a> <- cokrigagem(estima[[1]]$Estimativas, loc=gr,dados.comp=dados)</div><div><a href="http://preditos.gh">preditos.gh</a> <- volta.quad(<a href="http://md.cov.ck">md.cov.ck</a>,n.pontos=7,Variancia=FALSE)</div>
<div><a href="http://preditos.gh">preditos.gh</a> <- data.frame(<a href="http://preditos.gh">preditos.gh</a>)</div><div>write.table(<a href="http://preditos.gh">preditos.gh</a>,"pred_by4k7ns1000MBM.txt")</div>
<div>preditos.simu <- volta.cokri(<a href="http://md.cov.ck">md.cov.ck</a>,num.simu=1000,int.conf=0.95) #Ponto com problema quanto a processamento</div><div><br></div><div>preditos.simu <- data.frame(preditos.simu[[1]])</div>
<div>preditos.simu.ic <- data.frame(preditos.simu[[2]])</div><div>write.table(preditos.simu,"predsimu_by4k7ns1000MBM.txt")</div><div>write.table(preditos.simu.ic,"predsimuic_by4k7ns1000MBM.txt")</div>
<div>-------------------------------------------------------------------</div><div><br></div><div>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.</div>
<div><br></div><div>-------------------------------------------------------------------</div><div>volta.cokri <-</div><div>function(mat.cokri, num.simu, retorna.tudo=FALSE, int.conf=0.95){</div><div> nlinhas <- dim(mat.cokri[[1]])[1]</div>
<div> compos1 <- data.frame(matrix(nrow=nlinhas/2,ncol=3))</div><div> compos <- data.frame(matrix(nrow=nlinhas/2,ncol=3))</div><div><br></div><div> for(i in 1:num.simu){</div><div> g <- mvrnorm(n=1, mat.cokri[[1]],mat.cokri[[2]])</div>
<div> seq1 <- seq(1,nlinhas,by=2)</div><div> seq2 <- seq(2,nlinhas,by=2)</div><div> y1 <- g[seq1]</div><div> y2 <- g[seq2]</div><div> gerado <- data.frame(y1,y2)</div><div> compos <- agl(gerado)</div>
<div> compos1 <- cbind(compos,compos1)</div><div> }</div><div> </div><div> compos2 <- as.matrix(compos1)</div><div> </div><div> dim.vetor <- num.simu*3</div><div> sy1 <- seq(1,dim.vetor,by=3) </div>
<div> sy2 <- seq(2,dim.vetor,by=3) </div><div> sy3 <- seq(3,dim.vetor,by=3) </div><div> amostra1 <- as.matrix(compos2[,sy1],ncol=num.simu)</div><div> amostra2 <- as.matrix(compos2[,sy2],ncol=num.simu)</div>
<div> amostra3 <- as.matrix(compos2[,sy3],ncol=num.simu)</div><div><br></div><div>if(retorna.tudo == TRUE){</div><div> retorna <- list()</div><div> retorna[[1]] <- amostra1</div><div> retorna[[2]] <- amostra2</div>
<div> retorna[[3]] <- amostra3</div><div> return(retorna)</div><div>}</div><div><br></div><div>if(retorna.tudo == FALSE){</div><div> med1 <- apply(amostra1,1,mean)</div><div> med2 <- apply(amostra2,1,mean)</div>
<div> med3 <- apply(amostra3,1,mean)</div><div><br></div><div> q1 <- t(apply(amostra1,1,quantile,prob=c(1-int.conf,int.conf)))</div><div> q2 <- t(apply(amostra2,1,quantile,prob=c(1-int.conf,int.conf)))</div><div>
q3 <- t(apply(amostra3,1,quantile,prob=c(1-int.conf,int.conf)))</div><div> quantis <- cbind(q1,q2,q3)</div><div> resultado <- list()</div><div> resultado$preditos <- data.frame(med1,med2,med3)</div><div> names(resultado$preditos) <- c("Areia","Silte","Argila")</div>
<div> resultado$intervalo <- data.frame(quantis)</div><div> names(resultado$intervalo) <- c("LI Areia", "LS Areia", "LI Silte", "LS Silte", "LI Argila", "LS Argila")</div>
<div> return(resultado)}</div><div><br></div><div>}</div><div>-------------------------------------------------------------------</div><div><br></div><div>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.</div>
<div><br></div><div>A seguir está o trecho da função desmembrada.</div><div><br></div><div>-------------------------------------------------------------------</div><div>mat.cokri <- <a href="http://md.cov.ck">md.cov.ck</a> </div>
<div><br></div><div>retorna.tudo=FALSE</div><div>int.conf=0.95</div><div><br></div><div>nlinhas <- dim(mat.cokri[[1]])[1]</div><div>compos1 <- data.frame(matrix(nrow=nlinhas/2,ncol=3))</div><div>compos <- data.frame(matrix(nrow=nlinhas/2,ncol=3))</div>
<div><br></div><div>mvrnorm <- function(n = 1, mu, Sigma, tol=1e-6, empirical = FALSE)</div><div>{</div><div> p <- length(mu)</div><div> if(!all(dim(Sigma) == c(p,p))) stop("incompatible arguments")</div>
<div> eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE) </div><div> ev <- eS$values</div><div> if(!all(ev >= -tol*abs(ev[1L]))) ev <- ev *(-1)</div><div> X <- matrix(rnorm(p * n), n)</div><div>
if(empirical) {</div><div> X <- scale(X, TRUE, FALSE) </div><div> X <- X %*% svd(X, nu = 0)$v </div><div> X <- scale(X, FALSE, TRUE) </div><div> }</div><div> X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X)</div>
<div> nm <- names(mu)</div><div> if(is.null(nm) && !is.null(dn <- dimnames(Sigma))) nm <- dn[[1L]]</div><div> dimnames(X) <- list(nm, NULL)</div><div> if(n == 1) drop(X) else t(X)</div><div>
}</div><div><br></div><div> m1 <- mat.cokri[[2]][1:2601,1:2601] <span class="Apple-tab-span" style="white-space:pre"> </span>#divisao da matriz original</div><div> m2 <- mat.cokri[[2]][2602:5202,1:2601]</div>
<div> m3 <- mat.cokri[[2]][1:2601,2602:5202]</div><div> m4 <- mat.cokri[[2]][2602:5202,2602:5202]</div><div><br></div><div>#chamada a funcao mclapply (neste exemplo esta com 1:4, mas que na execução real seria 1:1000)</div>
<div>mclapply(1:4, function(i) {</div><div><br></div><div> g1 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m1)</div><div> g2 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m2)</div><div> g3 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m3)</div>
<div> g4 <- mvrnorm(n=1, mat.cokri[[1]][1:2601],m4)</div><div> list(g1, g2, g3, g4)</div><div><br></div><div>},mc.cores=4</div><div>)</div><div> g <- list(list(1:2601, g1), list(1:2601, g2), list(1:2601, g3),list(1:2601, g4))</div>
<div> seq1 <- seq(1,nlinhas,by=2)</div><div> seq2 <- seq(2,nlinhas,by=2)</div><div> y1 <- g[seq1]</div><div> y2 <- g[seq2]</div><div> gerado <- data.frame(y1,y2)</div><div> compos <- agl(gerado)</div>
<div> compos1 <- cbind(compos,compos1)</div><div><br></div><div> </div><div> compos2 <- as.matrix(compos1)</div><div> </div><div><br></div><div>-------------------------------------------------------------------</div>
<div><br></div><div>Esses trechos de código podem reproduzir o problema. Se não tiver o pacote geoComp, posso lhe enviar por email.</div><div>Meu email é <a href="mailto:beleti.junior@gmail.com">beleti.junior@gmail.com</a></div>
<div><br></div><div>Sobre a possivel solução apresentada, consegui aplica-lá.</div><div>Mais uma vez agradeço a sua disposição e peço desculpas pelo encomodo.</div><div><br></div><div>Att, Carlos Beleti.</div><br>