R help - Sep 2008 - Computation of contour values - Speeding up computation

Dear R useRs, 

i have the following code to compute values needed for a contour plot

############################################################

"myContour" <- function(a, b, plist, veca, vecb, dim)
{
  tmpb <- seq(0.5 * b, 1.5 * b, length=dim)
  tmpa  <- seq(0.5 * a, 1.5 * a, length=dim)

  z <- matrix(0, nrow=dim, ncol=dim)
  for(i in 1:dim)
  {
    for(j in 1:dim)
    {
      z[i, j] <- posteriorPdf(a=tmpa[j], b=tmpb[i], 
                              plist=plist, veca=veca, vecb=vecb) 
    }
  }
}

"posteriorPdf" <- function(a, b, plist, veca, vecb)
{
  res <- sum(plist[, 1] *       
             exp(vecb[, 1] * log(vecb[, 2]) + (vecb[, 1] - 1.0) * log(b) 
                 - vecb[, 2] * b - lgamma(vecb[, 1])) *
             exp(veca[, 1] * log(veca[, 2]) + (veca[, 1] - 1.0) * log(a) 
                 - veca[, 2] * a - lgamma(veca[, 1])))
          
  return(res)          
}

plist <- matrix(0, 100, 3)
plist[, 1] <- runif(100) 
veca <- vecb <- matrix(0, 100, 2)

veca[, 1] <- seq(20, 50, len=100)
veca[, 2] <- seq(10, 20, len=100)

vecb[, 1] <- seq(50, 200, len=100)
vecb[, 2] <- seq(1000, 400000, len=100)

myContour(a=20, b=0.01, plist=plist, veca=veca, vecb=vecb, dim=50)

############################################################

this is part of my other computations which i do with R. Here i recognized, that
my functions myContour and posteriorPdf took a long time of my computations. The
key to speed this up is to avoid the two for-loops in myContour, i think. I
tried a lot to do this with apply or something like that, but i didn't get
it.

If you have any advice how i can to this computations fast, i would be very
thankful, one idea is to use external c-code?


best regards
Andreas 
--

Andreas Wittmann wrote:
> Dear R useRs, 
> 
> i have the following code to compute values needed for a contour plot
> 
> ############################################################
> 
> "myContour" <- function(a, b, plist, veca, vecb, dim)
> {
>   tmpb <- seq(0.5 * b, 1.5 * b, length=dim)
>   tmpa  <- seq(0.5 * a, 1.5 * a, length=dim)
> 
>   z <- matrix(0, nrow=dim, ncol=dim)
>   for(i in 1:dim)
>   {
>     for(j in 1:dim)
>     {
>       z[i, j] <- posteriorPdf(a=tmpa[j], b=tmpb[i], 
>                               plist=plist, veca=veca, vecb=vecb) 
>     }
>   }
> }
> 
> "posteriorPdf" <- function(a, b, plist, veca, vecb)
> {
>   res <- sum(plist[, 1] *       
>              exp(vecb[, 1] * log(vecb[, 2]) + (vecb[, 1] - 1.0) * log(b) 
>                  - vecb[, 2] * b - lgamma(vecb[, 1])) *
>              exp(veca[, 1] * log(veca[, 2]) + (veca[, 1] - 1.0) * log(a) 
>                  - veca[, 2] * a - lgamma(veca[, 1])))
>           
>   return(res)          
> }
> 
> plist <- matrix(0, 100, 3)
> plist[, 1] <- runif(100) 
> veca <- vecb <- matrix(0, 100, 2)
> 
> veca[, 1] <- seq(20, 50, len=100)
> veca[, 2] <- seq(10, 20, len=100)
> 
> vecb[, 1] <- seq(50, 200, len=100)
> vecb[, 2] <- seq(1000, 400000, len=100)
> 
> myContour(a=20, b=0.01, plist=plist, veca=veca, vecb=vecb, dim=50)
> 
> ############################################################
> 
> this is part of my other computations which i do with R. Here i recognized,
that my functions myContour and posteriorPdf took a long time of my
computations. The key to speed this up is to avoid the two for-loops in
myContour, i think. I tried a lot to do this with apply or something like that,
but i didn't get it.
> 
> If you have any advice how i can to this computations fast, i would be very
thankful, one idea is to use external c-code?
> 

It takes 0.8 seconds on my machine. Not worth working on it, is it?

Your problem is that you are applying many calculations for all 
iterations of the inner loop, even if the result won't change, example:
lgamma(veca[, 1]) will be calculated dim^2 times!

Hence you can improve your loop considerably:



myContour <- function(a, b, plist, veca, vecb, dim)
{
   tmpb <- seq(0.5 * b, 1.5 * b, length=dim)
   tmpa  <- seq(0.5 * a, 1.5 * a, length=dim)

   z <- matrix(0, nrow=dim, ncol=dim)
   plist1 <- plist[, 1]
   vecb1l2 <- vecb[, 1] * log(vecb[, 2])
   vecb11 <- vecb[, 1] - 1
   vecb1lg <- lgamma(vecb[, 1])
   vecb2 <- vecb[, 2]
   veca1l2 <- veca[, 1] * log(veca[, 2])
   veca11 <- veca[, 1] - 1
   veca2 <- veca[, 2]
   veca1lg <- lgamma(veca[, 1])

   for(i in 1:dim)
   {
     for(j in 1:dim)
     {
       z[i, j] <- sum(plist1 * exp(vecb1l2 + vecb11 * log(tmpb[i]) - 
vecb2 * tmpb[i] - vecb1lg) * exp(veca1l2 + veca11 * log(tmpa[j]) - veca2 
* tmpa[j] - veca1lg))
     }
   }
   z
}

Uwe Ligges





> best regards
> Andreas 
> --
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.