The R Project SVN R

Rev

Rev 19233 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 19233 Rev 21689
Line 155... Line 155...
155 
	}
 155 
	}
156 
	UNPROTECT(1); /* x */
 156 
	UNPROTECT(1); /* x */
157 
    }
 157 
    }
158 
}
 158 
}
159 
 
 159 
 
- 
 
 160 
static void genptry(int n, double *p, double *ptry, double scale, void *ex)
- 
 
 161 
{
- 
 
 162 
    SEXP s, x;
- 
 
 163 
    int i;
- 
 
 164 
    OptStruct OS = (OptStruct) ex;
- 
 
 165 
    PROTECT_INDEX ipx;
- 
 
 166 
 
- 
 
 167 
    if (!isNull(OS->R_gcall)) {
- 
 
 168 
	/* user defined generation of candidate point */
- 
 
 169 
      	PROTECT(x = allocVector(REALSXP, n));
- 
 
 170 
	for (i = 0; i < n; i++) {
- 
 
 171 
	    if (!R_FINITE(p[i])) error("non-finite value supplied by optim");
- 
 
 172 
	    REAL(x)[i] = p[i] * (OS->parscale[i]);
- 
 
 173 
	}
- 
 
 174 
	SETCADR(OS->R_gcall, x);
- 
 
 175 
	PROTECT_WITH_INDEX(s = eval(OS->R_gcall, OS->R_env), &ipx);
- 
 
 176 
	REPROTECT(s = coerceVector(s, REALSXP), ipx);
- 
 
 177 
	if(LENGTH(s) != n)
- 
 
 178 
	    error("candidate point in optim evaluated to length %d not %d",
- 
 
 179 
		  LENGTH(s), n);
- 
 
 180 
	for (i = 0; i < n; i++)
- 
 
 181 
	    ptry[i] = REAL(s)[i] / (OS->parscale[i]);
- 
 
 182 
	UNPROTECT(2);
- 
 
 183 
    }
- 
 
 184 
    else {  /* default Gaussian Markov kernel */
- 
 
 185 
        for (i = 0; i < n; i++)
- 
 
 186 
            ptry[i] = p[i] + scale * norm_rand();  /* new candidate point */
- 
 
 187 
    }
- 
 
 188 
}
- 
 
 189 
 
160 
/* par fn gr method options */
 190 
/* par fn gr method options */
161 
SEXP do_optim(SEXP call, SEXP op, SEXP args, SEXP rho)
 191 
SEXP do_optim(SEXP call, SEXP op, SEXP args, SEXP rho)
162 
{
 192 
{
163 
    SEXP par, fn, gr, method, options, tmp, slower, supper;
 193 
    SEXP par, fn, gr, method, options, tmp, slower, supper;
164 
    SEXP res, value, counts, conv;
 194 
    SEXP res, value, counts, conv;
Line 219... Line 249...
219 
	for (i = 0; i < npar; i++)
 249 
	for (i = 0; i < npar; i++)
220 
	    REAL(par)[i] = opar[i] * (OS->parscale[i]);
 250 
	    REAL(par)[i] = opar[i] * (OS->parscale[i]);
221 
	grcount = NA_INTEGER;
 251 
	grcount = NA_INTEGER;
222 
 
 252 
 
223 
    }
 253 
    }
224 
    else if (strcmp(tn, "SANN") == 0) {
 254 
    else if (strcmp(tn, "SANN") == 0) {
225 
      tmax = asInteger(getListElement(options, "tmax"));
 255 
        tmax = asInteger(getListElement(options, "tmax"));
226 
      temp = asReal(getListElement(options, "temp"));
 256 
        temp = asReal(getListElement(options, "temp"));
227 
      if (tmax == NA_INTEGER) error("tmax is not an integer");
 257 
        if (tmax == NA_INTEGER) error("tmax is not an integer");
- 
 
 258 
        if (!isNull(gr)) {
- 
 
 259 
            if (!isFunction(gr)) error("gr is not a function");
- 
 
 260 
                PROTECT(OS->R_gcall = lang2(gr, R_NilValue));
- 
 
 261 
        } else {
- 
 
 262 
	    PROTECT(OS->R_gcall = R_NilValue); /* for balance */
- 
 
 263 
        }
228 
      samin (npar, dpar, &val, fminfn, maxit, tmax, temp, trace, (void *)OS);
 264 
        samin (npar, dpar, &val, fminfn, maxit, tmax, temp, trace, (void *)OS);
229 
      for (i = 0; i < npar; i++)
 265 
        for (i = 0; i < npar; i++)
230 
	  REAL(par)[i] = dpar[i] * (OS->parscale[i]);
 266 
            REAL(par)[i] = dpar[i] * (OS->parscale[i]);
231 
      fncount = maxit;
 267 
        fncount = maxit;
232 
      grcount = NA_INTEGER;
 268 
        grcount = NA_INTEGER;
- 
 
 269 
        UNPROTECT(1);  /* OS->R_gcall */
233 
 
 270 
 
234 
    } else if (strcmp(tn, "BFGS") == 0) {
 271 
    } else if (strcmp(tn, "BFGS") == 0) {
235 
	SEXP ndeps;
 272 
	SEXP ndeps;
236 
 
 273 
 
237 
	nREPORT = asInteger(getListElement(options, "REPORT"));
 274 
	nREPORT = asInteger(getListElement(options, "REPORT"));
Line 1038... Line 1075...
1038 
   evaluations at each temperature. Returned quantities are pb
 1075 
   evaluations at each temperature. Returned quantities are pb
1039 
   (the location of the minimum), and yb (the minimum value of
 1076 
   (the location of the minimum), and yb (the minimum value of
1040 
   the function func).  Author: Adrian Trapletti
 1077 
   the function func).  Author: Adrian Trapletti
1041 
*/
 1078 
*/
1042 
{
 1079 
{
1043 
    long i, j;
 1080 
    long j;
1044 
    int k, its, itdoc;
 1081 
    int k, its, itdoc;
1045 
    double t, y, dy, ytry, scale;
 1082 
    double t, y, dy, ytry, scale;
1046 
    double *p, *dp, *ptry;
 1083 
    double *p, *dp, *ptry;
1047 
 
 1084 
 
1048 
    p = vect (n); dp = vect (n); ptry = vect (n);
 1085 
    p = vect (n); dp = vect (n); ptry = vect (n);
Line 1056... Line 1093...
1056 
	Rprintf ("sann objective function values\n");
 1093 
	Rprintf ("sann objective function values\n");
1057 
	Rprintf ("initial       value %f\n", *yb);
 1094 
	Rprintf ("initial       value %f\n", *yb);
1058 
    }
 1095 
    }
1059 
    scale = 1.0/ti;
 1096 
    scale = 1.0/ti;
1060 
    its = itdoc = 1;
 1097 
    its = itdoc = 1;
1061 
    while (its < maxit) { /* cool down system */
 1098 
    while (its < maxit) {  /* cool down system */
1062 
	t = ti/log((double)its + E1);  /* temperature annealing schedule */
 1099 
	t = ti/log((double)its + E1);  /* temperature annealing schedule */
1063 
	k = 1;
 1100 
	k = 1;
1064 
	while ((k <= tmax) && (its < maxit))  /* iterate at constant temperature */
 1101 
	while ((k <= tmax) && (its < maxit))  /* iterate at constant temperature */
1065 
	{
 1102 
	{
1066 
	    for (i = 0; i < n; i++)
 - 
 
1067 
		dp[i] = scale * t * norm_rand();  /* random perturbation */
 - 
 
1068 
	    for (i = 0; i < n; i++)
 - 
 
1069 
		ptry[i] = p[i] + dp[i];  /* new candidate point */
 1103 
            genptry(n, p, ptry, scale * t, ex);  /* generate new candidate point */
1070 
	    ytry = fminfn (n, ptry, ex);
 1104 
	    ytry = fminfn (n, ptry, ex);
1071 
	    if (!R_FINITE(ytry)) ytry = big;
 1105 
	    if (!R_FINITE(ytry)) ytry = big;
1072 
	    dy = ytry - y;
 1106 
	    dy = ytry - y;
1073 
	    if ((dy <= 0.0) || (unif_rand() < exp(-dy/t))) {  /* accept new point? */
 1107 
	    if ((dy <= 0.0) || (unif_rand() < exp(-dy/t))) {  /* accept new point? */
1074 
		for (j = 0; j < n; j++) p[j] = ptry[j];
 1108 
		for (j = 0; j < n; j++) p[j] = ptry[j];