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