[Rd] API for optimization with Simulated annealing

[Jakub Dvorak]

Dear list,

I was trying to use the R API for optimization method "Simulated annealing"
    void samin(int n, double *x, double *Fmin, optimfn fn, int maxit,
               int tmax, double temp, int trace, void *ex);
but I encountered the following problem:

The implementation of the function samin (as seen in src/main/optim.c)
passes its void * argument "ex" into the function genptry (implemented
in the same source file) and the function genptry expects that the void *
argument "ex" points to a struct opt_struct, from which genptry takes a
function generating a new candidate point (if it is present). The struct
opt_struct seems to be intended only for internal use in implementation
of optim called from R. Using opt_struct in my C code calling samin
encounters the following problems:
1. The definition of opt_struct is present inside the file optim.c, thus
my source has no access to it.
2. The function genptry uses opt_struct attribute R_gcall of type SEXP which
doesn't allow to specify a function writen in C for generating a new
candidate point. This should be either an R function or NULL, meaning
that the default generating function is used.

Hence, a C user has to copy the definition of opt_struct to his
code and fill NULL to its R_gcall component.

This may cause problems, when opt_struct changes in future. Also,
it does not allow to use a user defined C function for generating
a new candidate point.

In order to solve the problems above, I would like to suggest that
samin does not pass the argument "ex" to genptry, which is fixed.
Instead, the function generating a new candidate point should be given
as an new pointer-to-function argument, say new_candidate, of function
samin. NULL value would mean "use the default". When called from R,
a (pointer to) function similar to genptry would be given as new_candidate
argument. Then, samin would call the function *new_candidate and pass
"ex" to it. Since "new_candidate" and "ex" are created in the same
environment, they may share the definition of an appropriate opt_struct.

Please, find a suggestion for a corresponding patch below.

The code modified according to my patch satifies the following general
requirement: The only thing that an implementation of the function
samin as an API should do with its argument void * ex is to pass it to
functions given by other arguments of samin. Then the user has full control
on usage of argument "ex" and is completely free to carry any data in it
because only he/she is responsible on right interpretation and use of the data.

Any change of API makes current users of the API unhappy. Still, I would
like to suggest the patch below for consideration, since, in my opinion, the
current samin API is not used frequently due to the problems described above.

        Jakub Dvorak

And here is the patch:

--- R-devel_2007-10-22/src/include/R_ext/Applic.h	2007-08-31 17:53:41.000000000 +0200
+++ R-devel_my_source/src/include/R_ext/Applic.h	2007-10-23 23:02:43.397529120 +0200
@@ -65,6 +65,7 @@
 /* main/optim.c */
 typedef double optimfn(int, double *, void *);
 typedef void optimgr(int, double *, double *, void *);
+typedef void sagenptry(int, double *, double *, double, void *);

 void vmmin(int n, double *b, double *Fmin,
 	   optimfn fn, optimgr gr, int maxit, int trace,
@@ -82,8 +83,8 @@
 	    double *Fmin, optimfn fn, optimgr gr, int *fail, void *ex,
 	    double factr, double pgtol, int *fncount, int *grcount,
 	    int maxit, char *msg, int trace, int nREPORT);
-void samin(int n, double *pb, double *yb, optimfn fn, int maxit,
-	   int tmax, double ti, int trace, void *ex);
+void samin(int n, double *pb, double *yb, optimfn fn, sagenptry *new_candidate,
+       int maxit, int tmax, double ti, int trace, void *ex);



--- R-devel_2007-10-22/src/main/optim.c	2007-08-15 17:50:12.000000000 +0200
+++ R-devel_my_source/src/main/optim.c	2007-10-23 23:02:43.397529120 +0200
@@ -171,6 +171,13 @@
     }
 }

+static void genptry_default(int n, double *p, double *ptry, double scale, void *ex)
+{  /* default Gaussian Markov kernel */
+    int i;
+    for (i = 0; i < n; i++)
+        ptry[i] = p[i] + scale * norm_rand();  /* new candidate point */
+}
+
 static void genptry(int n, double *p, double *ptry, double scale, void *ex)
 {
     SEXP s, x;
@@ -196,10 +203,8 @@
 	    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 */
-    }
+    else
+        genptry_default(n, p, ptry, scale, ex);
 }

 /* par fn gr method options */
@@ -275,7 +280,7 @@
         } else {
 	    PROTECT(OS->R_gcall = R_NilValue); /* for balance */
         }
-        samin (npar, dpar, &val, fminfn, maxit, tmax, temp, trace, (void *)OS);
+        samin (npar, dpar, &val, fminfn, &genptry, maxit, tmax, temp, trace, (void *)OS);
         for (i = 0; i < npar; i++)
             REAL(par)[i] = dpar[i] * (OS->parscale[i]);
         fncount = npar > 0 ? maxit : 1;
@@ -1084,8 +1089,8 @@
 #define E1 1.7182818  /* exp(1.0)-1.0 */
 #define STEPS 100

-void samin(int n, double *pb, double *yb, optimfn fminfn, int maxit,
-	   int tmax, double ti, int trace, void *ex)
+void samin(int n, double *pb, double *yb, optimfn fminfn, sagenptry *new_candidate,
+       int maxit, int tmax, double ti, int trace, void *ex)

 /* Given a starting point pb[0..n-1], simulated annealing minimization
    is performed on the function fminfn. The starting temperature
@@ -1121,9 +1126,11 @@
     while (its < maxit) {  /* cool down system */
 	t = ti/log((double)its + E1);  /* temperature annealing schedule */
 	k = 1;
+    if (new_candidate == NULL)
+        new_candidate = &genptry_default;
 	while ((k <= tmax) && (its < maxit))  /* iterate at constant temperature */
 	{
-            genptry(n, p, ptry, scale * t, ex);  /* generate new candidate point */
+        new_candidate(n, p, ptry, scale * t, ex);  /* generate new candidate point */
 	    ytry = fminfn (n, ptry, ex);
 	    if (!R_FINITE(ytry)) ytry = big;
 	    dy = ytry - y;