[R-br] Criação de função fortran em pacote (base)
Benilton Carvalho
beniltoncarvalho em gmail.com
Sexta Março 23 11:13:01 BRT 2012
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 em 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 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