From adrian.trapletti@lmttrading.com  Fri Oct 25 22:32:32 2002
From: adrian.trapletti@lmttrading.com (Adrian Trapletti)
Date: Fri, 25 Oct 2002 23:32:32 +0200
Subject: [R] Combinatorial Optimisation
References: <Pine.LNX.4.31.0210251259530.10820-100000@gannet.stats>
Message-ID: <3DB9B870.D7D565C1@lmttrading.com>

This is a multi-part message in MIME format.
--------------D871E36B191E8352EBC0EAA5
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

ripley@stats.ox.ac.uk wrote:

> Adrian,
>
> Thanks for the idea.  If you were willing to contribute a patch to do this
> we would be happy to incorporate it.

Patches for optim.c and optim.Rd against R-1.6.0 source are attached. I didn't test too extensively, but the basic things and the
traveling salesman example seem to work cleanly. I would be glad if you could review and maybe enhance, in particular the S related
code, in optim.c, because there I am really not an expert.

Furthermore, as I was scanning through optim.c I was not sure if line 99 (original)

     df[i] = REAL(s)[i] * (OS->parscale[i])/(OS->fnscale);

is correct. Shouldn't that be

    df[i] = REAL(s)[i] / ((OS->parscale[i])*(OS->fnscale));

to get back to the optim-internal scaling?

best
Adrian

--
Dr. Adrian Trapletti             Phone :             +41 (0) 1 994 5631
Trapletti Statistical Computing  Mobile:             +41 (0)76 370 5631
Wildsbergstrasse 31              Fax   :             +41 (0) 1 994 5631
CH-8610 Uster                    Email :  mailto:a.trapletti@bluewin.ch
Switzerland                      WWW   : http://trapletti.homelinux.com



--------------D871E36B191E8352EBC0EAA5
Content-Type: text/plain; charset=us-ascii;
 name="patch.optim.Rd"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="patch.optim.Rd"

7c7
<   box-constrained optimization.
---
>   box-constrained optimization and simulated annealing.
21,22c21,24
<    \code{"Nelder-Mead"} and \code{"SANN"} method. If it is \code{NULL}
<    and it is needed, a finite-difference approximation will be used.}
---
>    \code{"Nelder-Mead"} method. If it is \code{NULL} and it is needed, a
>    finite-difference approximation will be used. For the \code{"SANN"}
>    method it is a function to generate a new candidate point. If it is
>    \code{NULL} a default Gaussian Markov kernel is used.}
59c61
<   Method \code{"SANN"} is a variant of simulated annealing
---
>   Method \code{"SANN"} is by default a variant of simulated annealing
64,69c66,73
<   acceptance probability. The next candidate point is generated from a
<   Gaussian Markov kernel with scale proportional to the actual temperature.
<   Temperatures are decreased according to the logarithmic cooling
<   schedule as given in Belisle (1992, p. 890). Note that the
<   \code{"SANN"} method depends critically on the settings of the
<   control parameters.  It is not a general-purpose method but can be
---
>   acceptance probability. By default the next candidate point is
>   generated from a Gaussian Markov kernel with scale proportional to the
>   actual temperature. If a function to generate a new candidate point is
>   given, method \code{"SANN"} can also be used to solve combinatorial
>   optimization problems. Temperatures are decreased according to the
>   logarithmic cooling schedule as given in Belisle (1992, p. 890). Note
>   that the \code{"SANN"} method depends critically on the settings of
>   the control parameters. It is not a general-purpose method but can be 
176c180
< 
---
>   
240a245,290
> 
> ## Combinatorial optimization: Traveling salesman problem
> library(mva)
> 
> data(eurodist)
> eurodistmat <- as.matrix(eurodist)
> 
> distance <- function(sq) {  # Target function
>     sq2 <- embed(sq, 2)
>     return(sum(eurodistmat[cbind(sq2[,2],sq2[,1])]))
> }
> 
> genseq <- function(sq) {  # Generate new candidate sequence 
>     idx <- seq(2, NROW(eurodistmat)-1, by=1)
>     changepoints <- sample(idx, size=2, replace=FALSE)
>     tmp <- sq[changepoints[1]]
>     sq[changepoints[1]] <- sq[changepoints[2]]
>     sq[changepoints[2]] <- tmp
>     return(sq)
> }
> 
> sq <- c(1,2:NROW(eurodistmat),1)  # Initial sequence
> distance(sq)
> 
> res <- optim(sq, distance, genseq, method="SANN",
>              control=list(maxit=20000, temp=2000 ,trace=T))
> res  # Near optimum distance around 12842
> 
> loc <- cmdscale(eurodist)
> rx <- range(x <- loc[,1])
> ry <- range(y <- -loc[,2])
> tspinit <- loc[sq,]
> tspres <- loc[res$par,]
> s <- seq(NROW(tspres)-1)
> 
> plot(x, y, type="n", asp=1, xlab="", ylab="",
>      main="initial solution of traveling salesman problem")
> arrows(tspinit[s,1], -tspinit[s,2], tspinit[s+1,1], -tspinit[s+1,2],
>        angle=10, col="green") 
> text(x, y, names(eurodist), cex=0.8)
> 
> plot(x, y, type="n", asp=1, xlab="", ylab="",
>      main="optim() 'solving' traveling salesman problem") 
> arrows(tspres[s,1], -tspres[s,2], tspres[s+1,1], -tspres[s+1,2],
>        angle=10, col="red")
> text(x, y, names(eurodist), cex=0.8)

--------------D871E36B191E8352EBC0EAA5
Content-Type: text/plain; charset=us-ascii;
 name="patch.optim.c"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="patch.optim.c"

*** ./src/main/optim.c.origin	Mon Apr  8 08:19:48 2002
--- ./src/main/optim.c	Fri Oct 25 18:27:41 2002
***************
*** 155,164 ****
--- 155,193 ----
  	}
  	UNPROTECT(1); /* x */
      }
  }
  
+ static void genptry(int n, double *p, double *ptry, double scale, void *ex)
+ {    
+     SEXP s, x;
+     int i;
+     OptStruct OS = (OptStruct) ex;
+     PROTECT_INDEX ipx;
+ 
+     if (!isNull(OS->R_gcall)) {  /* user defined generation of candidate point */
+       	PROTECT(x = allocVector(REALSXP, n));
+ 	for (i = 0; i < n; i++) {
+ 	    if (!R_FINITE(p[i])) error("non-finite value supplied by optim");
+ 	    REAL(x)[i] = p[i] * (OS->parscale[i]);
+ 	}
+ 	SETCADR(OS->R_gcall, x);
+ 	PROTECT_WITH_INDEX(s = eval(OS->R_gcall, OS->R_env), &ipx);
+ 	REPROTECT(s = coerceVector(s, REALSXP), ipx);
+ 	if(LENGTH(s) != n)
+ 	    error("candidate point in optim evaluated to length %d not %d",
+ 		  LENGTH(s), n);
+ 	for (i = 0; i < n; i++)
+ 	    ptry[i] = REAL(s)[i] / (OS->parscale[i]);
+ 	UNPROTECT(2);
+     } 
+     else {  /* default Gaussian Markov kernel */
+         for (i = 0; i < n; i++)
+             ptry[i] = p[i] + scale * norm_rand();  /* new candidate point */
+     }
+ }
+ 
  /* par fn gr method options */
  SEXP do_optim(SEXP call, SEXP op, SEXP args, SEXP rho)
  {
      SEXP par, fn, gr, method, options, tmp, slower, supper;
      SEXP res, value, counts, conv;
***************
*** 219,237 ****
  	for (i = 0; i < npar; i++)
  	    REAL(par)[i] = opar[i] * (OS->parscale[i]);
  	grcount = NA_INTEGER;
  
      }
!     else if (strcmp(tn, "SANN") == 0) {
!       tmax = asInteger(getListElement(options, "tmax"));
!       temp = asReal(getListElement(options, "temp"));
!       if (tmax == NA_INTEGER) error("tmax is not an integer");
!       samin (npar, dpar, &val, fminfn, maxit, tmax, temp, trace, (void *)OS);
!       for (i = 0; i < npar; i++)
! 	  REAL(par)[i] = dpar[i] * (OS->parscale[i]);
!       fncount = maxit;
!       grcount = NA_INTEGER;
  
      } else if (strcmp(tn, "BFGS") == 0) {
  	SEXP ndeps;
  
  	nREPORT = asInteger(getListElement(options, "REPORT"));
--- 248,273 ----
  	for (i = 0; i < npar; i++)
  	    REAL(par)[i] = opar[i] * (OS->parscale[i]);
  	grcount = NA_INTEGER;
  
      }
!     else if (strcmp(tn, "SANN") == 0) {	
!         tmax = asInteger(getListElement(options, "tmax"));
!         temp = asReal(getListElement(options, "temp"));
!         if (tmax == NA_INTEGER) error("tmax is not an integer");
!         if (!isNull(gr)) {
!             if (!isFunction(gr)) error("gr is not a function");
!                 PROTECT(OS->R_gcall = lang2(gr, R_NilValue));
!         } else {
! 	    PROTECT(OS->R_gcall = R_NilValue); /* for balance */
!         }
!         samin (npar, dpar, &val, fminfn, maxit, tmax, temp, trace, (void *)OS);
!         for (i = 0; i < npar; i++)
!             REAL(par)[i] = dpar[i] * (OS->parscale[i]);
!         fncount = maxit;
!         grcount = NA_INTEGER;
!         UNPROTECT(1);  /* OS->R_gcall */
  
      } else if (strcmp(tn, "BFGS") == 0) {
  	SEXP ndeps;
  
  	nREPORT = asInteger(getListElement(options, "REPORT"));
***************
*** 1056,1074 ****
  	Rprintf ("sann objective function values\n");
  	Rprintf ("initial       value %f\n", *yb);
      }
      scale = 1.0/ti;
      its = itdoc = 1;
!     while (its < maxit) { /* cool down system */
  	t = ti/log((double)its + E1);  /* temperature annealing schedule */
  	k = 1;
  	while ((k <= tmax) && (its < maxit))  /* iterate at constant temperature */
  	{
! 	    for (i = 0; i < n; i++)
! 		dp[i] = scale * t * norm_rand();  /* random perturbation */
! 	    for (i = 0; i < n; i++)
! 		ptry[i] = p[i] + dp[i];  /* new candidate point */
  	    ytry = fminfn (n, ptry, ex);
  	    if (!R_FINITE(ytry)) ytry = big;
  	    dy = ytry - y;
  	    if ((dy <= 0.0) || (unif_rand() < exp(-dy/t))) {  /* accept new point? */
  		for (j = 0; j < n; j++) p[j] = ptry[j];
--- 1092,1107 ----
  	Rprintf ("sann objective function values\n");
  	Rprintf ("initial       value %f\n", *yb);
      }
      scale = 1.0/ti;
      its = itdoc = 1;
!     while (its < maxit) {  /* cool down system */
  	t = ti/log((double)its + E1);  /* temperature annealing schedule */
  	k = 1;
  	while ((k <= tmax) && (its < maxit))  /* iterate at constant temperature */
  	{
!             genptry(n, p, ptry, scale * t, ex);  /* generate new candidate point */
  	    ytry = fminfn (n, ptry, ex);
  	    if (!R_FINITE(ytry)) ytry = big;
  	    dy = ytry - y;
  	    if ((dy <= 0.0) || (unif_rand() < exp(-dy/t))) {  /* accept new point? */
  		for (j = 0; j < n; j++) p[j] = ptry[j];

--------------D871E36B191E8352EBC0EAA5--

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel 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-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


[R] Combinatorial Optimisation
Adrian Trapletti