Skip to content

Full enumeration, how can this be vectorized

10 messages · Daniel Hoppe, Uwe Ligges, Rafael A. Irizarry +3 more

Daniel Hoppe
# 
Hi all,

I am currently using the code snippet below to minimize a function which I
wasn't able to minimize properly with optim (cliffy surface, constraints,
poles). Obviously this iteration is probably the poorest way of implementing
such a function and it is expectedly slow. I've already written some
vectorized functions, but I'm afraid that I'm missing the "big picture" in
this case, I don't know how to start. Probably I could build a 3 x
(length(x) * length(y) * length(z)) matrix holding each possible parameter
combinations. With my current step-size and value range that would be
slightly more than 2,000,000 triples, so that should not be a problem. Then
comes the part I'm totally clueless about, that is applying the target
function to each of the triples in the matrix in an efficient way. From the
resulting vector I could then easily take the minimum value index and get
the parameters from the parameter matrix.

Could someone kindly give me a hint how this could be implemented?

Thanks a lot,

Daniel


Code Snippet:
fullEnumeration <- function(par, stepSize = c(.03, .03, .03))
{
   res <- Inf

   best.x <- -1
   best.y <- -1
   best.z <- -1

   stepx <- stepSize[1]
   stepy <- stepSize[2]
   stepz <- stepSize[3]
   x <- seq(.01, 2, by=stepx)
   y <- seq(.01, 2, by=stepy)
   z <- seq(.01, 15, by=stepz)

   for (i in 1:length(x))
    for (j in 1:length(y))
        for (k in 1:length(z))
        {
		# Pass the current best result to the function being minimized to allow
for early termination
            newRes <- fun(c(x[i], y[j], z[k]), par, res)
            if (!is.na(newRes) && newRes < res)
            {
                res = newRes
                best.x <- x[i]
                best.y <- y[j]
                best.z <- z[k]
            }

        }
   return (c(
        best.x,
        best.y,
        best.z))

}

fun <- function(x, par, currentBestValue=Inf)
{
  x = x[1]
  y = x[2]
  z = x[3]
  TMax <- length(par)

  result <- 0
  for (myT in 1:TMax)
  {
    result <- result +
    (
    ((x * z) / (y-1) * (1 - (z / (z + myT))^(y-1) ))
    - par[myT]
    )^2

    # Allow for short-cut evaluation if running in complete enumeration mode

    if (is.na(result) || currentBestValue < result)
    {
        return (Inf)
    }

  }
  return (result)

}



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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 ...

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?

Uwe Ligges
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Daniel Hoppe
# 
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.
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

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Uwe Ligges
#
Daniel Hoppe wrote:
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Rafael A. Irizarry
# 
hi! 

i need to submit electronic figures to a journal that satisfy the 
following requirement:

"Colour figures must be supplied in CMYK not RGB colours."


can i produce graphs in R that satisfy this?

thanks,
rafael


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Carlos Ortega
# 
Hi,

Well you can produce JPEG, PNG or BMP images directly from R and later you
can transform the image to the desired type of file format with another
second program (one that provide the feature you need is, ImageMagick).

Carlos.


-----Mensaje original-----
De: owner-r-help at stat.math.ethz.ch
[mailto:owner-r-help at stat.math.ethz.ch]En nombre de Rafael A. Irizarry
Enviado el: martes, 26 de noviembre de 2002 15:02
Para: r-help at stat.math.ethz.ch
Asunto: [R] CMYK


hi!

i need to submit electronic figures to a journal that satisfy the
following requirement:

"Colour figures must be supplied in CMYK not RGB colours."


can i produce graphs in R that satisfy this?

thanks,
rafael


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._


_____
The information in this email is confidential and it may not be
disclosed or used by anyone other than the addressee. If you are not the
intended recipient, any disclosure, copying, distribution or any action taken or
omitted is prohibited and may be unlawful.
Minorplanet cannot accept responsibility for the accuracy or completeness of
this email as it has been transmitted over a public network. If you suspect
that the email may have been amended, please call the sender.


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Brian Ripley
#
On Tue, 26 Nov 2002, Rafael A. Irizarry wrote:

            

      
      
      
    

        
Almost any way you like, then import into Adobe PhotoShop 6/7 and produce
CMYK output.

Given that `CMYK' is not a standard (you have to know what the numbers are
going to mean on the output device), the only reliable way I know is to
use a tool like PhotoShop that is set up to do colour profiling.  There
are all sorts of problems, like the restricted gamut of most CMYK devices.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Matej Cepl
    

    

      
      #
    

  


  

      
Carlos Ortega wrote:

            

      
      
      
    

        
On slightly related theme -- is there anyway how to make these
EPS pictures less bulky? I am writing a paper with couple of
pairs() and they are all 600k, 800k, etc. And it is really heck
of waiting to work with them (moreover, currently it is three
pages, but PDF file is already 1MB). However, I would prefer to
use EPS before PNG/BMP (which would be probably smaller and
faster) because of device independency.

Any thoughts on that?

	Matej

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Carlos Ortega
    

    

      
      #
    

  


  

      
Hi,

Yes, at least for JPEG, you can define the level of compression (see the
help file for that "jpeg") and the smaller it is the bigger is the file, but
the higher is the compression the lower is the quality of the picture. You
decide your satisfaction level. In any case, in terms of size, PNG's format
gets lower file sizes.

Carlos.

-----Mensaje original-----
De: owner-r-help at stat.math.ethz.ch
[mailto:owner-r-help at stat.math.ethz.ch]En nombre de Matej Cepl
Enviado el: martes, 26 de noviembre de 2002 16:45
Para: r-help at stat.math.ethz.ch
Asunto: Re: [R] CMYK

        
Carlos Ortega wrote:

            

      
      
      
    

        
On slightly related theme -- is there anyway how to make these
EPS pictures less bulky? I am writing a paper with couple of
pairs() and they are all 600k, 800k, etc. And it is really heck
of waiting to work with them (moreover, currently it is three
pages, but PDF file is already 1MB). However, I would prefer to
use EPS before PNG/BMP (which would be probably smaller and
faster) because of device independency.

Any thoughts on that?

	Matej

--
Matej Cepl, matej at ceplovi.cz,
Finger: 89EF 4BC6 288A BF43 1BAB  25C3 E09F EF25 D964 84AC
138 Highland Ave. #10, Somerville, Ma 02143, (617) 623-1488

The difference between death and taxes is death doesn't get worse
every time Congress meets
    -- Will Rogers

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._


_____
The information in this email is confidential and it may not be
disclosed or used by anyone other than the addressee. If you are not the
intended recipient, any disclosure, copying, distribution or any action taken or
omitted is prohibited and may be unlawful.
Minorplanet cannot accept responsibility for the accuracy or completeness of
this email as it has been transmitted over a public network. If you suspect
that the email may have been amended, please call the sender.


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Uwe Ligges
    

    

      
      #
    

  


  

      
Matej Cepl wrote:

            

      
      
      
    

        
These are vector format and each plotted data point, each line etc. is
represented in the file, so size can be modelled as lm(size ~
number.of.elements) (yes, including an intercept). ;-)

Uwe Ligges
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._