Re: [R-br] Criação de função fortran em pacote (base)
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.
geoComp nao esta' no CRAN... se, ao inves de mandar para o meu email, for possivel apontar para o link onde eu posso baixa-lo, melhor. Sobre o problema que vc descreve... Eu acho pouco provavel que o maior problema seja o mvrnorm()... Dois complicadores claros sao: 1) o for(i in 1:num.simu) e 2) o cbind(compos, compos1). A unica coisa q vc precisaria fazer seria reescrever a parte interna desse for() (ha' operacoes q nao precisam ser realizadas) e, depois, paralelizar o proprio for() (nao ha' pq mexer nas funcoes usadas dentro do for)... Apesar de que eu chego a arriscar que, uma vez reescrito a parte interna do for, o ganho em paralelizar o for() pode nao ser significativo. b 2012/3/6 Junior Beleti <beleti.junior@gmail.com>:
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.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Postando a msg do Carlos abaixo, afim de manter a discussao inteira dentro do topico inicialmente criado por ele proprio: Olá Benilton, segue link para download do pacote geoComp: www.din.uem.br/~pg45178/geocomp/geoComp.tar.gz Sobre o problema de desempenho, acredito que as complicações sejam no mvrnorm() pois é lá que o script leva mais tempo na execução. Continuo trabalhando para encontrar solução ou melhoria no código. Obrigado pela atenção. Carlos. 2012/3/6 Benilton Carvalho <beniltoncarvalho@gmail.com>:
geoComp nao esta' no CRAN... se, ao inves de mandar para o meu email, for possivel apontar para o link onde eu posso baixa-lo, melhor.
Sobre o problema que vc descreve... Eu acho pouco provavel que o maior problema seja o mvrnorm()... Dois complicadores claros sao: 1) o for(i in 1:num.simu) e 2) o cbind(compos, compos1).
A unica coisa q vc precisaria fazer seria reescrever a parte interna desse for() (ha' operacoes q nao precisam ser realizadas) e, depois, paralelizar o proprio for() (nao ha' pq mexer nas funcoes usadas dentro do for)... Apesar de que eu chego a arriscar que, uma vez reescrito a parte interna do for, o ganho em paralelizar o for() pode nao ser significativo.
b
2012/3/6 Junior Beleti <beleti.junior@gmail.com>:
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.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Carlos, eu nao estou em condicoes de tentar nada que me tome mais de 1 minuto por esses dias... Mas se sua impressao e' que o mvrnorm() eh o problema (o q ainda nao acho que seja), entao manda bala... Uma sugestao (isso vai depender do seu sistema operacional, que nao lembro se vc ja citou e, se sim, nao encontrei nas mensagens anteriores): instale uma versao do LAPACK/BLAS ja' com suporte a processadores com multiplos nucleos e, entao, "link" LAPACK/BLAS ao R... O mvrnorm, se bem me recordo, usa LAPACK para fazer sua magica... e se o LAPACK que vc estiver disponivel ja' tiver suporte paralelo, seus problemas estarao "resolvidos"... Nota: alterar o LAPACK/BLAS pode nao ser trivial, mas eh mais facil q vc reinventar a roda (leia-se mvrnorm). b 2012/3/23 Benilton Carvalho <beniltoncarvalho@gmail.com>:
Postando a msg do Carlos abaixo, afim de manter a discussao inteira dentro do topico inicialmente criado por ele proprio:
Olá Benilton,
segue link para download do pacote geoComp: www.din.uem.br/~pg45178/geocomp/geoComp.tar.gz
Sobre o problema de desempenho, acredito que as complicações sejam no mvrnorm() pois é lá que o script leva mais tempo na execução.
Continuo trabalhando para encontrar solução ou melhoria no código.
Obrigado pela atenção.
Carlos.
2012/3/6 Benilton Carvalho <beniltoncarvalho@gmail.com>:
geoComp nao esta' no CRAN... se, ao inves de mandar para o meu email, for possivel apontar para o link onde eu posso baixa-lo, melhor.
Sobre o problema que vc descreve... Eu acho pouco provavel que o maior problema seja o mvrnorm()... Dois complicadores claros sao: 1) o for(i in 1:num.simu) e 2) o cbind(compos, compos1).
A unica coisa q vc precisaria fazer seria reescrever a parte interna desse for() (ha' operacoes q nao precisam ser realizadas) e, depois, paralelizar o proprio for() (nao ha' pq mexer nas funcoes usadas dentro do for)... Apesar de que eu chego a arriscar que, uma vez reescrito a parte interna do for, o ganho em paralelizar o for() pode nao ser significativo.
b
2012/3/6 Junior Beleti <beleti.junior@gmail.com>:
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.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
participantes (2)
-
Benilton Carvalho -
Junior Beleti