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

[Junior Beleti]

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.
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20120306/c29da6c2/attachment.html>