AW: [R] Full enumeration, how can this be vectorized

[Uwe Ligges]

Daniel Hoppe wrote:
> 
> >
> > Without looking closer, I think this is a very nice example
> > for a piece
> > of code that should be implemented in either C or Fortran (on the one
> > hand it is easy to implement, on the other hand its "S interpreted"
> > execution takes quite a long time). OK, quite a few possible
> > improvements of the R code are obvious, but I don't think this is the
> > best point to start ...
> >
> 
> Probably you are right about C and Fortran. But coming from the combination
> of Java and Windows it's hard to overcome one's inhibitions and get started
> with C :-). Beside that I'm specifically trying to improve my R-skills a
> little because I wrote quite a lot of code for my thesis with R and
> therefore I'm curious for better "code patterns".  In my R-code I vectorized
> the function fun (thanks Patrick), seems to run somehow faster already.
> 
> > Nevertheless, the code looks like you are searching on a grid. Are you
> > really sure whether there is no method in optim() which
> > performs better?
> >
> 
> I tried them but was not successful yet. The search range is constrained,
> and Nelder-Mead doesn't use constraints. L-BFGS-B handles the constraints,
> but it isn't able to handle functions with poles. The other ones did not
> qualify for the same reasons.
> 
> I wonder if a grid-search is a very unusual case or if there might be a need
> for a fast and generic implementation of such an algorithm.
> 
> Daniel
> 


Just for fun two steps on the way to faster code:


fullEnumeration <- function(par, stepSize = c(.03, .03, .03))
{
   res <- Inf
   seqpar <- seq(along = par)
   best <- rep(-1, 3)
   x <- seq(.01, 2, by = stepSize[1])
   y <- seq(.01, 2, by = stepSize[2])
   z <- seq(.01, 15, by = stepSize[3])
   for (xi in x)
     for (yj in y)
       for (zk in z){
         # Pass the current best result to the function being minimized
to allow for early termination
         newRes <- fun(xi, yj, zk, par, seqpar, res)
         ## UL: I think (newRes < res) is more often FALSE than
!is.na(newRes), 
         ## UL: so turn around these two statements in that case:
         if (!is.na(newRes) && newRes < res){
                res <- newRes
                best <- c(xi, yj, zk)
         }
       }
   return (best)
}

fun <- function(x, y, z, par, seqpar, currentBestValue = Inf)
{ 
  ym <- y - 1
  result <- sum((((x*z) / ym * (1 - (z / (z + seqpar))^ym)) - par)^2)
  # Allow for short-cut evaluation if running in complete enumeration
mode
  if(is.na(result) || currentBestValue < result)
      return(Inf)
  return(result)
}

[Not tested!!!]
==========================================================

Next step: The following will be much faster, I guess [code not
tested!], but fun() is combined with fillEnumeration, which is not
reusable for arbitrary "funs" know....


fullEnumeration <- function(par, stepSize = c(.03, .03, .03))
{
   res <- Inf
   seqpar <- seq(along = par)
   best <- rep(-1, 3)
   x <- seq(.01, 2, by = stepSize[1])
   y <- seq(.01, 2, by = stepSize[2])
   z <- seq(.01, 15, by = stepSize[3])
   for (zk in z){
     zterm <- (zk / (zk + seqpar))
     for (yj in (y-1)){
       yterm <- yj / (1 - zterm^yj)
       for (xi in x){
         # Pass the current best result to the function being minimized
to allow for early termination
         newRes <- sum((((xi*zk) / yterm) - par)^2)
         ## UL: I think (newRes < res) is more often FALSE than
!is.na(newRes), 
         ## UL: so turn around these two statements in that case:
         if (!is.na(newRes) && newRes < res){
                res <- newRes
                best <- c(xi, yj, zk)
         }
       }
     }
   }
   return (best)
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._