[R-br] Criação de função fortran em pacote (base)
Benilton Carvalho
beniltoncarvalho em gmail.com
Sexta Março 23 11:04:05 BRT 2012
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 em 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 em 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 em 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 em 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.
Mais detalhes sobre a lista de discussão R-br