The R Project SVN R

Rev

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

Rev 45173 Rev 45446
Line 128... Line 128...
128 
		SETCADR(OS->R_fcall, x);
 128 
		SETCADR(OS->R_fcall, x);
129 
		REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);
 129 
		REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);
130 
		REPROTECT(s = coerceVector(s, REALSXP), ipx);
 130 
		REPROTECT(s = coerceVector(s, REALSXP), ipx);
131 
		val2 = REAL(s)[0]/(OS->fnscale);
 131 
		val2 = REAL(s)[0]/(OS->fnscale);
132 
		df[i] = (val1 - val2)/(2 * eps);
 132 
		df[i] = (val1 - val2)/(2 * eps);
133 
#define DO_df_x 							\
 133 
#define DO_df_x							\
134 
		if(!R_FINITE(df[i])) 					\
 134 
		if(!R_FINITE(df[i]))					\
135 
		    error(("non-finite finite-difference value [%d]"), i+1);\
 135 
		    error(("non-finite finite-difference value [%d]"), i+1);\
136 
		REAL(x)[i] = p[i] * (OS->parscale[i])
 136 
		REAL(x)[i] = p[i] * (OS->parscale[i])
137 
 
 137 
 
138 
		DO_df_x;
 138 
		DO_df_x;
139 
		UNPROTECT(1);
 139 
		UNPROTECT(1);
Line 178... Line 178...
178 
    OptStruct OS = (OptStruct) ex;
 178 
    OptStruct OS = (OptStruct) ex;
179 
    PROTECT_INDEX ipx;
 179 
    PROTECT_INDEX ipx;
180 
 
 180 
 
181 
    if (!isNull(OS->R_gcall)) {
 181 
    if (!isNull(OS->R_gcall)) {
182 
	/* user defined generation of candidate point */
 182 
	/* user defined generation of candidate point */
183 
      	PROTECT(x = allocVector(REALSXP, n));
 183 
	PROTECT(x = allocVector(REALSXP, n));
184 
	for (i = 0; i < n; i++) {
 184 
	for (i = 0; i < n; i++) {
185 
	    if (!R_FINITE(p[i]))
 185 
	    if (!R_FINITE(p[i]))
186 
		error(_("non-finite value supplied by optim"));
 186 
		error(_("non-finite value supplied by optim"));
187 
	    REAL(x)[i] = p[i] * (OS->parscale[i]);
 187 
	    REAL(x)[i] = p[i] * (OS->parscale[i]);
188 
	}
 188 
	}
Line 195... Line 195...
195 
	for (i = 0; i < n; i++)
 195 
	for (i = 0; i < n; i++)
196 
	    ptry[i] = REAL(s)[i] / (OS->parscale[i]);
 196 
	    ptry[i] = REAL(s)[i] / (OS->parscale[i]);
197 
	UNPROTECT(2);
 197 
	UNPROTECT(2);
198 
    }
 198 
    }
199 
    else {  /* default Gaussian Markov kernel */
 199 
    else {  /* default Gaussian Markov kernel */
200 
        for (i = 0; i < n; i++)
 200 
	for (i = 0; i < n; i++)
201 
            ptry[i] = p[i] + scale * norm_rand();  /* new candidate point */
 201 
	    ptry[i] = p[i] + scale * norm_rand();  /* new candidate point */
202 
    }
 202 
    }
203 
}
 203 
}
204 
 
 204 
 
205 
/* par fn gr method options */
 205 
/* par fn gr method options */
206 
SEXP attribute_hidden do_optim(SEXP call, SEXP op, SEXP args, SEXP rho)
 206 
SEXP attribute_hidden do_optim(SEXP call, SEXP op, SEXP args, SEXP rho)
Line 264... Line 264...
264 
	    REAL(par)[i] = opar[i] * (OS->parscale[i]);
 264 
	    REAL(par)[i] = opar[i] * (OS->parscale[i]);
265 
	grcount = NA_INTEGER;
 265 
	grcount = NA_INTEGER;
266 
 
 266 
 
267 
    }
 267 
    }
268 
    else if (strcmp(tn, "SANN") == 0) {
 268 
    else if (strcmp(tn, "SANN") == 0) {
269 
        tmax = asInteger(getListElement(options, "tmax"));
 269 
	tmax = asInteger(getListElement(options, "tmax"));
270 
        temp = asReal(getListElement(options, "temp"));
 270 
	temp = asReal(getListElement(options, "temp"));
271 
        if (trace) trace = asInteger(getListElement(options, "REPORT"));
 271 
	if (trace) trace = asInteger(getListElement(options, "REPORT"));
272 
        if (tmax == NA_INTEGER) error(_("'tmax' is not an integer"));
 272 
	if (tmax == NA_INTEGER) error(_("'tmax' is not an integer"));
273 
        if (!isNull(gr)) {
 273 
	if (!isNull(gr)) {
274 
            if (!isFunction(gr)) error(_("'gr' is not a function"));
 274 
	    if (!isFunction(gr)) error(_("'gr' is not a function"));
275 
                PROTECT(OS->R_gcall = lang2(gr, R_NilValue));
 275 
		PROTECT(OS->R_gcall = lang2(gr, R_NilValue));
276 
        } else {
 276 
	} else {
277 
	    PROTECT(OS->R_gcall = R_NilValue); /* for balance */
 277 
	    PROTECT(OS->R_gcall = R_NilValue); /* for balance */
278 
        }
 278 
	}
279 
        samin (npar, dpar, &val, fminfn, maxit, tmax, temp, trace, (void *)OS);
 279 
	samin (npar, dpar, &val, fminfn, maxit, tmax, temp, trace, (void *)OS);
280 
        for (i = 0; i < npar; i++)
 280 
	for (i = 0; i < npar; i++)
281 
            REAL(par)[i] = dpar[i] * (OS->parscale[i]);
 281 
	    REAL(par)[i] = dpar[i] * (OS->parscale[i]);
282 
        fncount = npar > 0 ? maxit : 1;
 282 
	fncount = npar > 0 ? maxit : 1;
283 
        grcount = NA_INTEGER;
 283 
	grcount = NA_INTEGER;
284 
        UNPROTECT(1);  /* OS->R_gcall */
 284 
	UNPROTECT(1);  /* OS->R_gcall */
285 
 
 285 
 
286 
    } else if (strcmp(tn, "BFGS") == 0) {
 286 
    } else if (strcmp(tn, "BFGS") == 0) {
287 
	SEXP ndeps;
 287 
	SEXP ndeps;
288 
 
 288 
 
289 
	nREPORT = asInteger(getListElement(options, "REPORT"));
 289 
	nREPORT = asInteger(getListElement(options, "REPORT"));
Line 1124... Line 1124...
1124 
    while (its < maxit) {  /* cool down system */
 1124 
    while (its < maxit) {  /* cool down system */
1125 
	t = ti/log((double)its + E1);  /* temperature annealing schedule */
 1125 
	t = ti/log((double)its + E1);  /* temperature annealing schedule */
1126 
	k = 1;
 1126 
	k = 1;
1127 
	while ((k <= tmax) && (its < maxit))  /* iterate at constant temperature */
 1127 
	while ((k <= tmax) && (its < maxit))  /* iterate at constant temperature */
1128 
	{
 1128 
	{
1129 
            genptry(n, p, ptry, scale * t, ex);  /* generate new candidate point */
 1129 
	    genptry(n, p, ptry, scale * t, ex);  /* generate new candidate point */
1130 
	    ytry = fminfn (n, ptry, ex);
 1130 
	    ytry = fminfn (n, ptry, ex);
1131 
	    if (!R_FINITE(ytry)) ytry = big;
 1131 
	    if (!R_FINITE(ytry)) ytry = big;
1132 
	    dy = ytry - y;
 1132 
	    dy = ytry - y;
1133 
	    if ((dy <= 0.0) || (unif_rand() < exp(-dy/t))) {  /* accept new point? */
 1133 
	    if ((dy <= 0.0) || (unif_rand() < exp(-dy/t))) {  /* accept new point? */
1134 
		for (j = 0; j < n; j++) p[j] = ptry[j];
 1134 
		for (j = 0; j < n; j++) p[j] = ptry[j];