The R Project SVN R

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 r 1 
/*
1160 maechler 2 
 *  R : A Computer Language for Statistical Data Analysis
2 r 3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
20764 maechler 4 
 *  Copyright (C) 1998--2002  Robert Gentleman, Ross Ihaka and the
5 
 *			      R Development Core Team
2 r 6 
 *
7 
 *  This program is free software; you can redistribute it and/or modify
8 
 *  it under the terms of the GNU General Public License as published by
9 
 *  the Free Software Foundation; either version 2 of the License, or
10 
 *  (at your option) any later version.
11 
 *
12 
 *  This program is distributed in the hope that it will be useful,
13 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 
 *  GNU General Public License for more details.
16 
 *
17 
 *  You should have received a copy of the GNU General Public License
18 
 *  along with this program; if not, write to the Free Software
5458 ripley 19 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
2 r 20 
 */
21 
 
5187 hornik 22 
#ifdef HAVE_CONFIG_H
7701 hornik 23 
#include <config.h>
5187 hornik 24 
#endif
25 
 
11499 ripley 26 
#include <Defn.h>
27 
#include <Print.h>		/* for printRealVector() */
28 
#include <Rmath.h>
29 
#include <R_ext/Applic.h>
30 
#include <R_ext/RS.h>			/* for Memcpy */
2 r 31 
 
32 
 
25717 maechler 33 
/* One Dimensional Minimization --- just wrapper for
34 
 * Brent's "fmin" --> ../appl/fmin.c */
2 r 35 
 
7762 ripley 36 
struct callinfo {
37 
  SEXP R_fcall;
38 
  SEXP R_env;
39 
} ;
2 r 40 
 
7762 ripley 41 
/*static SEXP R_fcall1;
42 
  static SEXP R_env1; */
43 
 
44 
static double fcn1(double x, struct callinfo *info)
2 r 45 
{
1881 ihaka 46 
    SEXP s;
7762 ripley 47 
    REAL(CADR(info->R_fcall))[0] = x;
48 
    s = eval(info->R_fcall, info->R_env);
1881 ihaka 49 
    switch(TYPEOF(s)) {
50 
    case INTSXP:
51 
	if (length(s) != 1) goto badvalue;
52 
	if (INTEGER(s)[0] == NA_INTEGER) {
4165 rgentlem 53 
	    warning("NA replaced by maximum positive value");
1881 ihaka 54 
	    return DBL_MAX;
2 r 55 
	}
1881 ihaka 56 
	else return INTEGER(s)[0];
57 
	break;
58 
    case REALSXP:
59 
	if (length(s) != 1) goto badvalue;
5107 maechler 60 
	if (!R_FINITE(REAL(s)[0])) {
4165 rgentlem 61 
	    warning("NA/Inf replaced by maximum positive value");
1881 ihaka 62 
	    return DBL_MAX;
63 
	}
64 
	else return REAL(s)[0];
65 
	break;
66 
    default:
67 
	goto badvalue;
68 
    }
69 
 badvalue:
7762 ripley 70 
    error("invalid function value in 'optimize'");
1881 ihaka 71 
    return 0;/* for -Wall */
2 r 72 
}
73 
 
2610 maechler 74 
/* fmin(f, xmin, xmax tol) */
2 r 75 
SEXP do_fmin(SEXP call, SEXP op, SEXP args, SEXP rho)
76 
{
1881 ihaka 77 
    double xmin, xmax, tol;
7762 ripley 78 
    SEXP v, res;
79 
    struct callinfo info;
2 r 80 
 
1881 ihaka 81 
    checkArity(op, args);
82 
    PrintDefaults(rho);
2 r 83 
 
1881 ihaka 84 
    /* the function to be minimized */
2 r 85 
 
1881 ihaka 86 
    v = CAR(args);
87 
    if (!isFunction(v))
5731 ripley 88 
	errorcall(call, "attempt to minimize non-function");
1881 ihaka 89 
    args = CDR(args);
2 r 90 
 
1881 ihaka 91 
    /* xmin */
2 r 92 
 
1881 ihaka 93 
    xmin = asReal(CAR(args));
5107 maechler 94 
    if (!R_FINITE(xmin))
6191 maechler 95 
	errorcall(call, "invalid xmin value");
1881 ihaka 96 
    args = CDR(args);
2 r 97 
 
1881 ihaka 98 
    /* xmax */
2 r 99 
 
1881 ihaka 100 
    xmax = asReal(CAR(args));
5107 maechler 101 
    if (!R_FINITE(xmax))
6191 maechler 102 
	errorcall(call, "invalid xmax value");
1881 ihaka 103 
    if (xmin >= xmax)
5731 ripley 104 
	errorcall(call, "xmin not less than xmax");
1881 ihaka 105 
    args = CDR(args);
2 r 106 
 
1881 ihaka 107 
    /* tol */
2 r 108 
 
1881 ihaka 109 
    tol = asReal(CAR(args));
5107 maechler 110 
    if (!R_FINITE(tol) || tol <= 0.0)
6191 maechler 111 
	errorcall(call, "invalid tol value");
2 r 112 
 
7762 ripley 113 
    info.R_env = rho;
114 
    PROTECT(info.R_fcall = lang2(v, R_NilValue));
115 
    PROTECT(res = allocVector(REALSXP, 1));
10172 luke 116 
    SETCADR(info.R_fcall, allocVector(REALSXP, 1));
20764 maechler 117 
    REAL(res)[0] = Brent_fmin(xmin, xmax,
7762 ripley 118 
			      (double (*)(double, void*)) fcn1, &info, tol);
119 
    UNPROTECT(2);
120 
    return res;
2 r 121 
}
122 
 
123 
 
124 
 
2610 maechler 125 
/* One Dimensional Root Finding --  just wrapper code for Brent's "zeroin" */
2 r 126 
 
127 
 
7762 ripley 128 
static double fcn2(double x, struct callinfo *info)
2 r 129 
{
1881 ihaka 130 
    SEXP s;
7762 ripley 131 
    REAL(CADR(info->R_fcall))[0] = x;
132 
    s = eval(info->R_fcall, info->R_env);
1881 ihaka 133 
    switch(TYPEOF(s)) {
134 
    case INTSXP:
135 
	if (length(s) != 1) goto badvalue;
136 
	if (INTEGER(s)[0] == NA_INTEGER) {
4165 rgentlem 137 
	    warning("NA replaced by maximum positive value");
1881 ihaka 138 
	    return	DBL_MAX;
2 r 139 
	}
1881 ihaka 140 
	else return INTEGER(s)[0];
141 
	break;
142 
    case REALSXP:
143 
	if (length(s) != 1) goto badvalue;
5107 maechler 144 
	if (!R_FINITE(REAL(s)[0])) {
4165 rgentlem 145 
	    warning("NA/Inf replaced by maximum positive value");
1881 ihaka 146 
	    return DBL_MAX;
147 
	}
148 
	else return REAL(s)[0];
149 
	break;
150 
    default:
151 
	goto badvalue;
152 
    }
153 
 badvalue:
5731 ripley 154 
    error("invalid function value in 'zeroin'");
1881 ihaka 155 
    return 0;/* for -Wall */
1016 maechler 156 
 
2 r 157 
}
158 
 
2612 maechler 159 
/* zeroin(f, xmin, xmax, tol, maxiter) */
2 r 160 
SEXP do_zeroin(SEXP call, SEXP op, SEXP args, SEXP rho)
161 
{
1881 ihaka 162 
    double xmin, xmax, tol;
2612 maechler 163 
    int iter;
164 
    SEXP v, res;
4275 pd 165 
    struct callinfo info;
2 r 166 
 
1881 ihaka 167 
    checkArity(op, args);
168 
    PrintDefaults(rho);
2 r 169 
 
1881 ihaka 170 
    /* the function to be minimized */
2 r 171 
 
1881 ihaka 172 
    v = CAR(args);
173 
    if (!isFunction(v))
5731 ripley 174 
	errorcall(call,"attempt to minimize non-function");
1881 ihaka 175 
    args = CDR(args);
2 r 176 
 
1881 ihaka 177 
    /* xmin */
2 r 178 
 
1881 ihaka 179 
    xmin = asReal(CAR(args));
5107 maechler 180 
    if (!R_FINITE(xmin))
6191 maechler 181 
	errorcall(call, "invalid xmin value");
1881 ihaka 182 
    args = CDR(args);
2 r 183 
 
1881 ihaka 184 
    /* xmax */
2 r 185 
 
1881 ihaka 186 
    xmax = asReal(CAR(args));
5107 maechler 187 
    if (!R_FINITE(xmax))
6191 maechler 188 
	errorcall(call, "invalid xmax value");
1881 ihaka 189 
    if (xmin >= xmax)
5731 ripley 190 
	errorcall(call, "xmin not less than xmax");
1881 ihaka 191 
    args = CDR(args);
2 r 192 
 
1881 ihaka 193 
    /* tol */
2 r 194 
 
1881 ihaka 195 
    tol = asReal(CAR(args));
5107 maechler 196 
    if (!R_FINITE(tol) || tol <= 0.0)
6191 maechler 197 
	errorcall(call, "invalid tol value");
2612 maechler 198 
    args = CDR(args);
2 r 199 
 
2612 maechler 200 
    /* maxiter */
201 
    iter = asInteger(CAR(args));
202 
    if (iter <= 0)
5731 ripley 203 
	errorcall(call, "maxiter must be positive");
2612 maechler 204 
 
7762 ripley 205 
    info.R_env = rho;
206 
    PROTECT(info.R_fcall = lang2(v, R_NilValue)); /* the info used in fcn2() */
10172 luke 207 
    SETCADR(info.R_fcall, allocVector(REALSXP, 1));
2612 maechler 208 
    PROTECT(res = allocVector(REALSXP, 3));
209 
    REAL(res)[0] =
7758 ripley 210 
	R_zeroin(xmin, xmax,   (double (*)(double, void*)) fcn2,
211 
		 (void *) &info, &tol, &iter);
2612 maechler 212 
    REAL(res)[1] = (double)iter;
213 
    REAL(res)[2] = tol;
214 
    UNPROTECT(2);
215 
    return res;
2 r 216 
}
3813 maechler 217 
 
2 r 218 
 
219 
 
1881 ihaka 220 
/* General Nonlinear Optimization */
2 r 221 
 
7166 bates 222 
#define FT_SIZE 5		/* default size of table to store computed
223 
				   function values */
2 r 224 
 
7166 bates 225 
typedef struct {
226 
  double   fval;
227 
  double  *x;
228 
  double  *grad;
229 
  double  *hess;
230 
} ftable;
2 r 231 
 
7166 bates 232 
typedef struct {
20764 maechler 233 
  SEXP R_fcall;	      /* unevaluated call to R function */
234 
  SEXP R_env;	      /* where to evaluate the calls */
7166 bates 235 
  int have_gradient;
236 
  int have_hessian;
237 
/*  int n;	      -* length of the parameter (x) vector */
20764 maechler 238 
  int FT_size;	      /* size of table to store computed
7166 bates 239 
			 function values */
20764 maechler 240 
  int FT_last;	      /* Newest entry in the table */
7166 bates 241 
  ftable *Ftable;
242 
} function_info;
6455 bates 243 
 
244 
/* Initialize the storage in the table of computed function values */
245 
 
7166 bates 246 
static void FT_init(int n, int FT_size, function_info *state)
6455 bates 247 
{
248 
    int i, j;
7166 bates 249 
    int have_gradient, have_hessian;
250 
    ftable *Ftable;
6623 maechler 251 
 
7166 bates 252 
    have_gradient = state->have_gradient;
253 
    have_hessian = state->have_hessian;
254 
 
255 
    Ftable = (ftable *)R_alloc(FT_size, sizeof(ftable));
256 
 
257 
    for (i = 0; i < FT_size; i++) {
258 
	Ftable[i].x = (double *)R_alloc(n, sizeof(double));
6455 bates 259 
				/* initialize to unlikely parameter values */
260 
	for (j = 0; j < n; j++) {
261 
	    Ftable[i].x[j] = DBL_MAX;
262 
	}
263 
	if (have_gradient) {
7166 bates 264 
	    Ftable[i].grad = (double *)R_alloc(n, sizeof(double));
6455 bates 265 
	    if (have_hessian) {
7176 bates 266 
		Ftable[i].hess = (double *)R_alloc(n * n, sizeof(double));
6455 bates 267 
	    }
268 
	}
269 
    }
7166 bates 270 
    state->Ftable = Ftable;
271 
    state->FT_size = FT_size;
272 
    state->FT_last = -1;
6455 bates 273 
}
274 
 
275 
/* Store an entry in the table of computed function values */
276 
 
7176 bates 277 
static void FT_store(int n, const double f, const double *x, const double *grad,
7166 bates 278 
		     const double *hess, function_info *state)
6455 bates 279 
{
6623 maechler 280 
    int ind;
6455 bates 281 
 
7176 bates 282 
    ind = (++(state->FT_last)) % (state->FT_size);
7166 bates 283 
    state->Ftable[ind].fval = f;
284 
    Memcpy(state->Ftable[ind].x, x, n);
6455 bates 285 
    if (grad) {
20764 maechler 286 
	Memcpy(state->Ftable[ind].grad, grad, n);
6455 bates 287 
	if (hess) {
7166 bates 288 
	    Memcpy(state->Ftable[ind].hess, hess, n * n);
6455 bates 289 
	}
290 
    }
291 
}
292 
 
293 
/* Check for stored values in the table of computed function values.
294 
   Returns the index in the table or -1 for failure */
295 
 
7166 bates 296 
static int FT_lookup(int n, const double *x, function_info *state)
6455 bates 297 
{
298 
    double *ftx;
299 
    int i, j, ind, matched;
7166 bates 300 
    int FT_size, FT_last;
301 
    ftable *Ftable;
20764 maechler 302 
 
7166 bates 303 
    FT_last = state->FT_last;
304 
    FT_size = state->FT_size;
305 
    Ftable = state->Ftable;
6455 bates 306 
 
7166 bates 307 
    for (i = 0; i < FT_size; i++) {
308 
	ind = (FT_last - i) % FT_size;
6455 bates 309 
				/* why can't they define modulus correctly */
7166 bates 310 
	if (ind < 0) ind += FT_size;
6455 bates 311 
	ftx = Ftable[ind].x;
312 
	if (ftx) {
313 
	    matched = 1;
314 
	    for (j = 0; j < n; j++) {
20764 maechler 315 
		if (x[j] != ftx[j]) {
316 
		    matched = 0;
6455 bates 317 
		    break;
20764 maechler 318 
		}
6455 bates 319 
	    }
320 
	    if (matched) return ind;
20764 maechler 321 
	}
6455 bates 322 
    }
323 
    return -1;
324 
}
325 
 
1881 ihaka 326 
/* This how the optimizer sees them */
2 r 327 
 
7166 bates 328 
static void fcn(int n, const double x[], double *f, function_info
329 
		*state)
2 r 330 
{
7166 bates 331 
    SEXP s, R_fcall;
332 
    ftable *Ftable;
6455 bates 333 
    double *g = (double *) 0, *h = (double *) 0;
1881 ihaka 334 
    int i;
6623 maechler 335 
 
7166 bates 336 
    R_fcall = state->R_fcall;
337 
    Ftable = state->Ftable;
338 
    if ((i = FT_lookup(n, x, state)) >= 0) {
6455 bates 339 
	*f = Ftable[i].fval;
7166 bates 340 
	return;
6455 bates 341 
    }
342 
				/* calculate for a new value of x */
7166 bates 343 
    s = CADR(R_fcall);
8316 ripley 344 
    for (i = 0; i < n; i++) {
345 
	if (!R_FINITE(x[i])) error("non-finite value supplied by nlm");
1881 ihaka 346 
	REAL(s)[i] = x[i];
8316 ripley 347 
    }
7166 bates 348 
    s = eval(state->R_fcall, state->R_env);
1881 ihaka 349 
    switch(TYPEOF(s)) {
350 
    case INTSXP:
351 
	if (length(s) != 1) goto badvalue;
352 
	if (INTEGER(s)[0] == NA_INTEGER) {
4165 rgentlem 353 
	    warning("NA replaced by maximum positive value");
1881 ihaka 354 
	    *f = DBL_MAX;
2 r 355 
	}
1881 ihaka 356 
	else *f = INTEGER(s)[0];
357 
	break;
358 
    case REALSXP:
359 
	if (length(s) != 1) goto badvalue;
5107 maechler 360 
	if (!R_FINITE(REAL(s)[0])) {
4165 rgentlem 361 
	    warning("NA/Inf replaced by maximum positive value");
1881 ihaka 362 
	    *f = DBL_MAX;
363 
	}
364 
	else *f = REAL(s)[0];
365 
	break;
366 
    default:
367 
	goto badvalue;
368 
    }
7166 bates 369 
    if (state->have_gradient) {
370 
	g = REAL(coerceVector(getAttrib(s, install("gradient")), REALSXP));
371 
	if (state->have_hessian) {
372 
	    h = REAL(coerceVector(getAttrib(s, install("hessian")), REALSXP));
6455 bates 373 
	}
374 
    }
7166 bates 375 
    FT_store(n, *f, x, g, h, state);
376 
    return;
377 
 
1881 ihaka 378 
 badvalue:
5731 ripley 379 
    error("invalid function value in 'nlm' optimizer");
2 r 380 
}
381 
 
382 
 
7166 bates 383 
static void Cd1fcn(int n, const double x[], double *g, function_info *state)
2 r 384 
{
6455 bates 385 
    int ind;
386 
 
7166 bates 387 
    if ((ind = FT_lookup(n, x, state)) < 0) {	/* shouldn't happen */
388 
	fcn(n, x, g, state);
389 
	if ((ind = FT_lookup(n, x, state)) < 0) {
6455 bates 390 
	    error("function value caching for optimization is seriously confused.\n");
391 
	}
392 
    }
7166 bates 393 
    Memcpy(g, state->Ftable[ind].grad, n);
2 r 394 
}
395 
 
396 
 
7166 bates 397 
static void Cd2fcn(int nr, int n, const double x[], double *h,
398 
		   function_info *state)
2 r 399 
{
6455 bates 400 
    int j, ind;
401 
 
7166 bates 402 
    if ((ind = FT_lookup(n, x, state)) < 0) {	/* shouldn't happen */
403 
	fcn(n, x, h, state);
404 
	if ((ind = FT_lookup(n, x, state)) < 0) {
6455 bates 405 
	    error("function value caching for optimization is seriously confused.\n");
406 
	}
407 
    }
7166 bates 408 
    for (j = 0; j < n; j++) {  /* fill in lower triangle only */
20764 maechler 409 
	Memcpy( h + j*(n + 1), state->Ftable[ind].hess + j*(n + 1), n - j);
6455 bates 410 
    }
2 r 411 
}
412 
 
413 
 
414 
static double *fixparam(SEXP p, int *n, SEXP call)
415 
{
1881 ihaka 416 
    double *x;
417 
    int i;
2 r 418 
 
1881 ihaka 419 
    if (!isNumeric(p))
5731 ripley 420 
	errorcall(call, "numeric parameter expected");
2 r 421 
 
1881 ihaka 422 
    if (*n) {
423 
	if (LENGTH(p) != *n)
5731 ripley 424 
	    errorcall(call, "conflicting parameter lengths");
1881 ihaka 425 
    }
426 
    else {
427 
	if (LENGTH(p) <= 0)
5731 ripley 428 
	    errorcall(call, "invalid parameter length");
1881 ihaka 429 
	*n = LENGTH(p);
430 
    }
431 
 
432 
    x = (double*)R_alloc(*n, sizeof(double));
433 
    switch(TYPEOF(p)) {
434 
    case LGLSXP:
435 
    case INTSXP:
436 
	for (i = 0; i < *n; i++) {
437 
	    if (INTEGER(p)[i] == NA_INTEGER)
5731 ripley 438 
		errorcall(call, "missing value in parameter");
1881 ihaka 439 
	    x[i] = INTEGER(p)[i];
2 r 440 
	}
1881 ihaka 441 
	break;
442 
    case REALSXP:
443 
	for (i = 0; i < *n; i++) {
5107 maechler 444 
	    if (!R_FINITE(REAL(p)[i]))
6191 maechler 445 
		errorcall(call, "missing value in parameter");
1881 ihaka 446 
	    x[i] = REAL(p)[i];
2 r 447 
	}
1881 ihaka 448 
	break;
449 
    default:
5731 ripley 450 
	errorcall(call, "invalid parameter type");
1881 ihaka 451 
    }
452 
    return x;
2 r 453 
}
454 
 
455 
 
456 
static void invalid_na(SEXP call)
457 
{
5731 ripley 458 
    errorcall(call, "invalid NA value in parameter");
2 r 459 
}
460 
 
461 
 
462 
	/* Fatal errors - we don't deliver an answer */
463 
 
1016 maechler 464 
static void opterror(int nerr)
2 r 465 
{
1881 ihaka 466 
    switch(nerr) {
467 
    case -1:
5731 ripley 468 
	error("non-positive number of parameters in nlm");
1881 ihaka 469 
    case -2:
5731 ripley 470 
	error("nlm is inefficient for 1-d problems");
1881 ihaka 471 
    case -3:
5731 ripley 472 
	error("illegal gradient tolerance in nlm");
1881 ihaka 473 
    case -4:
5731 ripley 474 
	error("illegal iteration limit in nlm");
1881 ihaka 475 
    case -5:
5731 ripley 476 
	error("minimization function has no good digits in nlm");
1881 ihaka 477 
    case -6:
5731 ripley 478 
	error("no analytic gradient to check in nlm!");
1881 ihaka 479 
    case -7:
5731 ripley 480 
	error("no analytic Hessian to check in nlm!");
1881 ihaka 481 
    case -21:
5731 ripley 482 
	error("probable coding error in analytic gradient");
1881 ihaka 483 
    case -22:
5731 ripley 484 
	error("probable coding error in analytic Hessian");
2612 maechler 485 
    default:
5731 ripley 486 
	error("*** unknown error message (msg = %d) in nlm()\n*** should not happen!", nerr);
1881 ihaka 487 
    }
2 r 488 
}
489 
 
490 
 
491 
	/* Warnings - we return a value, but print a warning */
492 
 
523 maechler 493 
static void optcode(int code)
2 r 494 
{
1881 ihaka 495 
    switch(code) {
496 
    case 1:
497 
	Rprintf("Relative gradient close to zero.\n");
498 
	Rprintf("Current iterate is probably solution.\n");
499 
	break;
500 
    case 2:
501 
	Rprintf("Successive iterates within tolerance.\n");
502 
	Rprintf("Current iterate is probably solution.\n");
503 
	break;
504 
    case 3:
505 
	Rprintf("Last global step failed to locate a point lower than x.\n");
506 
	Rprintf("Either x is an approximate local minimum of the function,\n");
507 
	Rprintf("the function is too non-linear for this algorithm,\n");
3805 ripley 508 
	Rprintf("or steptol is too large.\n");
1881 ihaka 509 
	break;
510 
    case 4:
511 
	Rprintf("Iteration limit exceeded.  Algorithm failed.\n");
512 
	break;
513 
    case 5:
514 
	Rprintf("Maximum step size exceeded 5 consecutive times.\n");
515 
	Rprintf("Either the function is unbounded below,\n");
516 
	Rprintf("becomes asymptotic to a finite value\n");
517 
	Rprintf("from above in some direction,\n");
518 
	Rprintf("or stepmx is too small.\n");
519 
	break;
520 
    }
521 
    Rprintf("\n");
2 r 522 
}
523 
 
20764 maechler 524 
/* NOTE: The actual Dennis-Schnabel algorithm `optif9' is in ../appl/uncmin.c */
525 
 
2 r 526 
SEXP do_nlm(SEXP call, SEXP op, SEXP args, SEXP rho)
527 
{
7166 bates 528 
    SEXP value, names, v, R_gradientSymbol, R_hessianSymbol;
2 r 529 
 
1881 ihaka 530 
    double *x, *typsiz, fscale, gradtl, stepmx,
3813 maechler 531 
	steptol, *xpls, *gpls, fpls, *a, *wrk, dlt;
2 r 532 
 
7185 bates 533 
    int code, i, j, k, itnlim, method, iexp, omsg, msg,
1881 ihaka 534 
	n, ndigit, iagflg, iahflg, want_hessian, itncnt;
2 r 535 
 
1881 ihaka 536 
    char *vmax;
2 r 537 
 
20764 maechler 538 
/* .Internal(
539 
 *	nlm(function(x) f(x, ...), p, hessian, typsize, fscale,
540 
 *	    msg, ndigit, gradtol, stepmax, steptol, iterlim)
541 
 */
7166 bates 542 
    function_info *state;
543 
 
1881 ihaka 544 
    checkArity(op, args);
545 
    PrintDefaults(rho);
546 
    vmax = vmaxget();
2 r 547 
 
7166 bates 548 
    state = (function_info *) R_alloc(1, sizeof(function_info));
549 
 
1881 ihaka 550 
    /* the function to be minimized */
2 r 551 
 
7166 bates 552 
    state->R_env = rho;
1881 ihaka 553 
    v = CAR(args);
554 
    if (!isFunction(v))
5731 ripley 555 
	error("attempt to minimize non-function");
7166 bates 556 
    PROTECT(state->R_fcall = lang2(v, R_NilValue));
1881 ihaka 557 
    args = CDR(args);
2 r 558 
 
20764 maechler 559 
    /* `p' : inital parameter value */
2 r 560 
 
1881 ihaka 561 
    n = 0;
562 
    x = fixparam(CAR(args), &n, call);
563 
    args = CDR(args);
2 r 564 
 
20764 maechler 565 
    /* `hessian' : H. required? */
523 maechler 566 
 
1881 ihaka 567 
    want_hessian = asLogical(CAR(args));
568 
    if (want_hessian == NA_LOGICAL) want_hessian = 0;
569 
    args = CDR(args);
2 r 570 
 
20764 maechler 571 
    /* `typsize' : typical size of parameter elements */
2 r 572 
 
1881 ihaka 573 
    typsiz = fixparam(CAR(args), &n, call);
574 
    args = CDR(args);
2 r 575 
 
20764 maechler 576 
    /* `fscale' : expected function size */
2 r 577 
 
1881 ihaka 578 
    fscale = asReal(CAR(args));
5107 maechler 579 
    if (ISNA(fscale)) invalid_na(call);
1881 ihaka 580 
    args = CDR(args);
2 r 581 
 
20764 maechler 582 
    /* `msg' (bit pattern) */
1881 ihaka 583 
    omsg = msg = asInteger(CAR(args));
584 
    if (msg == NA_INTEGER) invalid_na(call);
585 
    args = CDR(args);
2 r 586 
 
1881 ihaka 587 
    ndigit = asInteger(CAR(args));
588 
    if (ndigit == NA_INTEGER) invalid_na(call);
589 
    args = CDR(args);
2 r 590 
 
1881 ihaka 591 
    gradtl = asReal(CAR(args));
5107 maechler 592 
    if (ISNA(gradtl)) invalid_na(call);
1881 ihaka 593 
    args = CDR(args);
2 r 594 
 
1881 ihaka 595 
    stepmx = asReal(CAR(args));
5107 maechler 596 
    if (ISNA(stepmx)) invalid_na(call);
1881 ihaka 597 
    args = CDR(args);
2 r 598 
 
3813 maechler 599 
    steptol = asReal(CAR(args));
5107 maechler 600 
    if (ISNA(steptol)) invalid_na(call);
1881 ihaka 601 
    args = CDR(args);
2 r 602 
 
20764 maechler 603 
    /* `iterlim' (def. 100) */
1881 ihaka 604 
    itnlim = asInteger(CAR(args));
605 
    if (itnlim == NA_INTEGER) invalid_na(call);
606 
    args = CDR(args);
2 r 607 
 
6455 bates 608 
    /* force one evaluation to check for the gradient and hessian */
609 
    iagflg = 0;			/* No analytic gradient */
610 
    iahflg = 0;			/* No analytic hessian */
7166 bates 611 
    state->have_gradient = 0;
612 
    state->have_hessian = 0;
6455 bates 613 
    R_gradientSymbol = install("gradient");
614 
    R_hessianSymbol = install("hessian");
615 
 
616 
    v = allocVector(REALSXP, n);
617 
    for (i = 0; i < n; i++) {
618 
	REAL(v)[i] = x[i];
619 
    }
10172 luke 620 
    SETCADR(state->R_fcall, v);
7166 bates 621 
    value = eval(state->R_fcall, state->R_env);
6455 bates 622 
 
623 
    v = getAttrib(value, R_gradientSymbol);
11347 ripley 624 
    if (v != R_NilValue) {
625 
	if (LENGTH(v) == n && (isReal(v) || isInteger(v))) {
626 
	    iagflg = 1;
627 
	    state->have_gradient = 1;
628 
	    v = getAttrib(value, R_hessianSymbol);
20764 maechler 629 
 
11347 ripley 630 
	    if (v != R_NilValue) {
631 
		if (LENGTH(v) == (n * n) && (isReal(v) || isInteger(v))) {
632 
		    iahflg = 1;
633 
		    state->have_hessian = 1;
634 
		} else {
635 
		    warning("hessian supplied is of the wrong length or mode, so ignored");
636 
		}
637 
	    }
638 
	} else {
639 
	    warning("gradient supplied is of the wrong length or mode, so ignored");
640 
	}
6455 bates 641 
    }
7176 bates 642 
    if (((msg/4) % 2) && !iahflg) { /* skip check of analytic Hessian */
6615 bates 643 
      msg -= 4;
644 
    }
7176 bates 645 
    if (((msg/2) % 2) && !iagflg) { /* skip check of analytic gradient */
6615 bates 646 
      msg -= 2;
647 
    }
7166 bates 648 
    FT_init(n, FT_SIZE, state);
1881 ihaka 649 
    /* Plug in the call to the optimizer here */
2 r 650 
 
1881 ihaka 651 
    method = 1;	/* Line Search */
7166 bates 652 
    iexp = iahflg ? 0 : 1; /* Function calls are expensive */
1881 ihaka 653 
    dlt = 1.0;
2 r 654 
 
1881 ihaka 655 
    xpls = (double*)R_alloc(n, sizeof(double));
656 
    gpls = (double*)R_alloc(n, sizeof(double));
657 
    a = (double*)R_alloc(n*n, sizeof(double));
658 
    wrk = (double*)R_alloc(8*n, sizeof(double));
2 r 659 
 
1881 ihaka 660 
    /*
2605 maechler 661 
     *	 Dennis + Schnabel Minimizer
1881 ihaka 662 
     *
2605 maechler 663 
     *	  SUBROUTINE OPTIF9(NR,N,X,FCN,D1FCN,D2FCN,TYPSIZ,FSCALE,
664 
     *	 +	   METHOD,IEXP,MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR,
3813 maechler 665 
     *	 +	   DLT,GRADTL,STEPMX,STEPTOL,
2605 maechler 666 
     *	 +	   XPLS,FPLS,GPLS,ITRMCD,A,WRK)
1881 ihaka 667 
     *
668 
     *
2605 maechler 669 
     *	 Note: I have figured out what msg does.
670 
     *	 It is actually a sum of bit flags as follows
671 
     *	   1 = don't check/warn for 1-d problems
672 
     *	   2 = don't check analytic gradients
673 
     *	   4 = don't check analytic hessians
674 
     *	   8 = don't print start and end info
675 
     *	  16 = print at every iteration
676 
     *	 Using msg=9 is absolutely minimal
677 
     *	 I think we always check gradients and hessians
1881 ihaka 678 
     */
2 r 679 
 
7166 bates 680 
    optif9(n, n, x, (fcn_p) fcn, (fcn_p) Cd1fcn, (d2fcn_p) Cd2fcn,
681 
	   state, typsiz, fscale, method, iexp, &msg, ndigit, itnlim,
682 
	   iagflg, iahflg, dlt, gradtl, stepmx, steptol, xpls, &fpls,
683 
	   gpls, &code, a, wrk, &itncnt);
2 r 684 
 
1881 ihaka 685 
    if (msg < 0)
686 
	opterror(msg);
687 
    if (code != 0 && (omsg&8) == 0)
688 
	optcode(code);
2 r 689 
 
1881 ihaka 690 
    if (want_hessian) {
691 
	PROTECT(value = allocVector(VECSXP, 6));
692 
	PROTECT(names = allocVector(STRSXP, 6));
7166 bates 693 
	fdhess(n, xpls, fpls, (fcn_p) fcn, state, a, n, &wrk[0], &wrk[n],
694 
	       ndigit, typsiz);
1881 ihaka 695 
	for (i = 0; i < n; i++)
696 
	    for (j = 0; j < i; j++)
697 
		a[i + j * n] = a[j + i * n];
698 
    }
699 
    else {
700 
	PROTECT(value = allocVector(VECSXP, 5));
701 
	PROTECT(names = allocVector(STRSXP, 5));
702 
    }
703 
    k = 0;
2 r 704 
 
10172 luke 705 
    SET_STRING_ELT(names, k, mkChar("minimum"));
706 
    SET_VECTOR_ELT(value, k, allocVector(REALSXP, 1));
707 
    REAL(VECTOR_ELT(value, k))[0] = fpls;
1881 ihaka 708 
    k++;
2 r 709 
 
10172 luke 710 
    SET_STRING_ELT(names, k, mkChar("estimate"));
711 
    SET_VECTOR_ELT(value, k, allocVector(REALSXP, n));
1881 ihaka 712 
    for (i = 0; i < n; i++)
10172 luke 713 
	REAL(VECTOR_ELT(value, k))[i] = xpls[i];
1881 ihaka 714 
    k++;
2 r 715 
 
10172 luke 716 
    SET_STRING_ELT(names, k, mkChar("gradient"));
717 
    SET_VECTOR_ELT(value, k, allocVector(REALSXP, n));
1881 ihaka 718 
    for (i = 0; i < n; i++)
10172 luke 719 
	REAL(VECTOR_ELT(value, k))[i] = gpls[i];
1881 ihaka 720 
    k++;
2 r 721 
 
1881 ihaka 722 
    if (want_hessian) {
10172 luke 723 
	SET_STRING_ELT(names, k, mkChar("hessian"));
724 
	SET_VECTOR_ELT(value, k, allocMatrix(REALSXP, n, n));
1881 ihaka 725 
	for (i = 0; i < n * n; i++)
10172 luke 726 
	    REAL(VECTOR_ELT(value, k))[i] = a[i];
1881 ihaka 727 
	k++;
728 
    }
2 r 729 
 
10172 luke 730 
    SET_STRING_ELT(names, k, mkChar("code"));
731 
    SET_VECTOR_ELT(value, k, allocVector(INTSXP, 1));
732 
    INTEGER(VECTOR_ELT(value, k))[0] = code;
1881 ihaka 733 
    k++;
523 maechler 734 
 
1881 ihaka 735 
    /* added by Jim K Lindsey */
10172 luke 736 
    SET_STRING_ELT(names, k, mkChar("iterations"));
737 
    SET_VECTOR_ELT(value, k, allocVector(INTSXP, 1));
738 
    INTEGER(VECTOR_ELT(value, k))[0] = itncnt;
1881 ihaka 739 
    k++;
740 
 
741 
    setAttrib(value, R_NamesSymbol, names);
742 
    vmaxset(vmax);
743 
    UNPROTECT(3);
744 
    return value;
2 r 745 
}