[R] Combinatorial Optimisation
Adrian Trapletti
adrian.trapletti at lmttrading.com
Fri Oct 25 23:32:32 CEST 2002
ripley at 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 at bluewin.ch
Switzerland WWW : http://trapletti.homelinux.com
-------------- next part --------------
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)
-------------- next part --------------
*** ./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];
More information about the R-help
mailing list