[R-br] Criação de função fortran em pacote (base)

[Benilton Carvalho]

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.