The R Project SVN – R – Diff – /trunk/src/main/optim.c


/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1999-2002  the R Development Core Team
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */
 
#include <Defn.h>
#include <R_ext/Random.h>	/* for the random number generation in
				   samin() */
#include <R_ext/Applic.h>	/* setulb() */
 
static SEXP getListElement(SEXP list, char *str)
{
    SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);
    int i;
 
    for (i = 0; i < length(list); i++)
	if (strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {
	    elmt = VECTOR_ELT(list, i);
	    break;
	}
    return elmt;
}
 
static double * vect(int n)
{
    return (double *)R_alloc(n, sizeof(double));
}
 
typedef struct opt_struct
{
    SEXP R_fcall;    /* function */
    SEXP R_gcall;    /* gradient */
    SEXP R_env;      /* where to evaluate the calls */
    double* ndeps;   /* tolerances for numerical derivatives */
    double fnscale;  /* scaling for objective */
    double* parscale;/* scaling for parameters */
    int usebounds;
    double* lower, *upper;
} opt_struct, *OptStruct;
 
 
 
static double fminfn(int n, double *p, void *ex)
{
    SEXP s, x;
    int i;
    double val;
    OptStruct OS = (OptStruct) ex;
    PROTECT_INDEX ipx;
 
    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_fcall, x);
    PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);
    REPROTECT(s = coerceVector(s, REALSXP), ipx);
    val = REAL(s)[0]/(OS->fnscale);
    UNPROTECT(2);
    return val;
}
 
static void fmingr(int n, double *p, double *df, void *ex)
{
    SEXP s, x;
    int i;
    double val1, val2, eps, epsused, tmp;
    OptStruct OS = (OptStruct) ex;
    PROTECT_INDEX ipx;
 
    if (!isNull(OS->R_gcall)) { /* analytical derivatives */
	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("gradient in optim evaluated to length %d not %d",
		  LENGTH(s), n);
	for (i = 0; i < n; i++)
	    df[i] = REAL(s)[i] * (OS->parscale[i])/(OS->fnscale);
	UNPROTECT(2);
    } else { /* numerical derivatives */
	PROTECT(x = allocVector(REALSXP, n));
	for (i = 0; i < n; i++) REAL(x)[i] = p[i] * (OS->parscale[i]);
	SETCADR(OS->R_fcall, x);
	if(OS->usebounds == 0) {
	    for (i = 0; i < n; i++) {
		eps = OS->ndeps[i];
		REAL(x)[i] = (p[i] + eps) * (OS->parscale[i]);
		SETCADR(OS->R_fcall, x);
		PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);
		REPROTECT(s = coerceVector(s, REALSXP), ipx);
		val1 = REAL(s)[0]/(OS->fnscale);
		REAL(x)[i] = (p[i] - eps) * (OS->parscale[i]);
		SETCADR(OS->R_fcall, x);
		REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);
		REPROTECT(s = coerceVector(s, REALSXP), ipx);
		val2 = REAL(s)[0]/(OS->fnscale);
		df[i] = (val1 - val2)/(2 * eps);
#define DO_df_x 							\
		if(!R_FINITE(df[i])) 					\
		    error("non-finite finite-difference value [%d]", i);\
		REAL(x)[i] = p[i] * (OS->parscale[i])
 
		DO_df_x;
		UNPROTECT(1);
	    }
	} else { /* usebounds */
	    for (i = 0; i < n; i++) {
		epsused = eps = OS->ndeps[i];
		tmp = p[i] + eps;
		if (tmp > OS->upper[i]) {
		    tmp = OS->upper[i];
		    epsused = tmp - p[i] ;
		}
		REAL(x)[i] = tmp * (OS->parscale[i]);
		SETCADR(OS->R_fcall, x);
		PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);
		REPROTECT(s = coerceVector(s, REALSXP), ipx);
		val1 = REAL(s)[0]/(OS->fnscale);
		tmp = p[i] - eps;
		if (tmp < OS->lower[i]) {
		    tmp = OS->lower[i];
		    eps = p[i] - tmp;
		}
		REAL(x)[i] = tmp * (OS->parscale[i]);
		SETCADR(OS->R_fcall, x);
		REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);
		REPROTECT(s = coerceVector(s, REALSXP), ipx);
		val2 = REAL(s)[0]/(OS->fnscale);
		df[i] = (val1 - val2)/(epsused + eps);
 
		DO_df_x;
		UNPROTECT(1);
	    }
	}
	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;
    int i, npar=0, *mask, trace, maxit, fncount, grcount, nREPORT, tmax;
    int ifail = 0;
    double *dpar, *opar, val, abstol, reltol, temp;
    char *tn;
    OptStruct OS;
    char *vmax;
 
    checkArity(op, args);
    vmax = vmaxget();
    OS = (OptStruct) R_alloc(1, sizeof(opt_struct));
    OS->usebounds = 0;
    OS->R_env = rho;
    par = CAR(args);
    args = CDR(args); fn = CAR(args);
    if (!isFunction(fn)) errorcall(call, "fn is not a function");
    args = CDR(args); gr = CAR(args);
    args = CDR(args); method = CAR(args);
    if (!isString(method)|| LENGTH(method) != 1)
	errorcall(call, "invalid method argument");
    tn = CHAR(STRING_ELT(method, 0));
    args = CDR(args); options = CAR(args);
    PROTECT(OS->R_fcall = lang2(fn, R_NilValue));
    PROTECT(par = coerceVector(duplicate(par), REALSXP));
    npar = LENGTH(par);
    dpar = vect(npar);
    opar = vect(npar);
    trace = asInteger(getListElement(options, "trace"));
    OS->fnscale = asReal(getListElement(options, "fnscale"));
    tmp = getListElement(options, "parscale");
    if (LENGTH(tmp) != npar)
	errorcall(call, "parscale is of the wrong length");
    PROTECT(tmp = coerceVector(tmp, REALSXP));
    OS->parscale = vect(npar);
    for (i = 0; i < npar; i++) OS->parscale[i] = REAL(tmp)[i];
    UNPROTECT(1);
    for (i = 0; i < npar; i++)
	dpar[i] = REAL(par)[i] / (OS->parscale[i]);
    PROTECT(res = allocVector(VECSXP, 5));
    PROTECT(value = allocVector(REALSXP, 1));
    PROTECT(counts = allocVector(INTSXP, 2));
    PROTECT(conv = allocVector(INTSXP, 1));
    abstol = asReal(getListElement(options, "abstol"));
    reltol = asReal(getListElement(options, "reltol"));
    maxit = asInteger(getListElement(options, "maxit"));
    if (maxit == NA_INTEGER) error("maxit is not an integer");
 
    if (strcmp(tn, "Nelder-Mead") == 0) {
	double alpha, beta, gamm;
 
	alpha = asReal(getListElement(options, "alpha"));
	beta = asReal(getListElement(options, "beta"));
	gamm = asReal(getListElement(options, "gamma"));
	nmmin(npar, dpar, opar, &val, fminfn, &ifail, abstol, reltol,
	      (void *)OS, alpha, beta, gamm, trace, &fncount, maxit);
	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"));
	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 */
	    ndeps = getListElement(options, "ndeps");
	    if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");
	    OS->ndeps = vect(npar);
	    PROTECT(ndeps = coerceVector(ndeps, REALSXP));
	    for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];
	    UNPROTECT(1);
	}
	mask = (int *) R_alloc(npar, sizeof(int));
	for (i = 0; i < npar; i++) mask[i] = 1;
	vmmin(npar, dpar, &val, fminfn, fmingr, maxit, trace, mask, abstol,
	      reltol, nREPORT, (void *)OS, &fncount, &grcount, &ifail);
	for (i = 0; i < npar; i++)
	    REAL(par)[i] = dpar[i] * (OS->parscale[i]);
	UNPROTECT(1); /* OS->R_gcall */
    } else if (strcmp(tn, "CG") == 0) {
	int type;
	SEXP ndeps;
 
	type = asInteger(getListElement(options, "type"));
	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 */
	    ndeps = getListElement(options, "ndeps");
	    if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");
	    OS->ndeps = vect(npar);
	    PROTECT(ndeps = coerceVector(ndeps, REALSXP));
	    for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];
	    UNPROTECT(1);
	}
	cgmin(npar, dpar, opar, &val, fminfn, fmingr, &ifail, abstol,
	      reltol, (void *)OS, type, trace, &fncount, &grcount, maxit);
	for (i = 0; i < npar; i++)
	    REAL(par)[i] = opar[i] * (OS->parscale[i]);
	UNPROTECT(1); /* OS->R_gcall */
 
    } else if (strcmp(tn, "L-BFGS-B") == 0) {
	SEXP ndeps, smsg;
	double *lower = vect(npar), *upper = vect(npar);
	int lmm, *nbd = (int *) R_alloc(npar, sizeof(int));
	double factr, pgtol;
	char msg[60];
 
	nREPORT = asInteger(getListElement(options, "REPORT"));
	factr = asReal(getListElement(options, "factr"));
	pgtol = asReal(getListElement(options, "pgtol"));
	lmm = asInteger(getListElement(options, "lmm"));
	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 */
	    ndeps = getListElement(options, "ndeps");
	    if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");
	    OS->ndeps = vect(npar);
	    PROTECT(ndeps = coerceVector(ndeps, REALSXP));
	    for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];
	    UNPROTECT(1);
	}
	args = CDR(args); slower = CAR(args); /* coerce in calling code */
	args = CDR(args); supper = CAR(args);
	for (i = 0; i < npar; i++) {
	    lower[i] = REAL(slower)[i] / (OS->parscale[i]);
	    upper[i] = REAL(supper)[i] / (OS->parscale[i]);
	    if (!R_FINITE(lower[i])) {
		if (!R_FINITE(upper[i])) nbd[i] = 0; else nbd[i] = 3;
	    } else {
		if (!R_FINITE(upper[i])) nbd[i] = 1; else nbd[i] = 2;
	    }
	}
	OS->usebounds = 1;
	OS->lower = lower;
	OS->upper = upper;
	lbfgsb(npar, lmm, dpar, lower, upper, nbd, &val, fminfn, fmingr,
	       &ifail, (void *)OS, factr, pgtol, &fncount, &grcount,
	       maxit, msg, trace, nREPORT);
	for (i = 0; i < npar; i++)
	    REAL(par)[i] = dpar[i] * (OS->parscale[i]);
	UNPROTECT(1); /* OS->R_gcall */
	PROTECT(smsg = allocVector(STRSXP, 1));
	SET_STRING_ELT(smsg, 0, mkChar(msg));
	SET_VECTOR_ELT(res, 4, smsg);
	UNPROTECT(1);
    } else
	errorcall(call, "unknown method");
 
    REAL(value)[0] = val * (OS->fnscale);
    SET_VECTOR_ELT(res, 0, par); SET_VECTOR_ELT(res, 1, value);
    INTEGER(counts)[0] = fncount; INTEGER(counts)[1] = grcount;
    SET_VECTOR_ELT(res, 2, counts);
    INTEGER(conv)[0] = ifail;
    SET_VECTOR_ELT(res, 3, conv);
    vmaxset(vmax);
    UNPROTECT(6);
    return res;
}
 
/* par fn gr options */
SEXP do_optimhess(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    SEXP par, fn, gr, options, tmp, ndeps, ans;
    OptStruct OS;
    int npar, i , j;
    double *dpar, *df1, *df2, eps;
    char *vmax;
 
    checkArity(op, args);
    vmax = vmaxget();
    OS = (OptStruct) R_alloc(1, sizeof(opt_struct));
    OS->usebounds = 0;
    OS->R_env = rho;
    par = CAR(args);
    npar = LENGTH(par);
    args = CDR(args); fn = CAR(args);
    if (!isFunction(fn)) errorcall(call, "fn is not a function");
    args = CDR(args); gr = CAR(args);
    args = CDR(args); options = CAR(args);
    OS->fnscale = asReal(getListElement(options, "fnscale"));
    tmp = getListElement(options, "parscale");
    if (LENGTH(tmp) != npar)
	errorcall(call, "parscale is of the wrong length");
    PROTECT(tmp = coerceVector(tmp, REALSXP));
    OS->parscale = vect(npar);
    for (i = 0; i < npar; i++) OS->parscale[i] = REAL(tmp)[i];
    UNPROTECT(1);
    PROTECT(OS->R_fcall = lang2(fn, R_NilValue));
    PROTECT(par = coerceVector(par, REALSXP));
    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 */
    }
    ndeps = getListElement(options, "ndeps");
    if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");
    OS->ndeps = vect(npar);
    PROTECT(ndeps = coerceVector(ndeps, REALSXP));
    for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];
    UNPROTECT(1);
    PROTECT(ans = allocMatrix(REALSXP, npar, npar));
    dpar = vect(npar);
    for (i = 0; i < npar; i++)
	dpar[i] = REAL(par)[i] / (OS->parscale[i]);
    df1 = vect(npar);
    df2 = vect(npar);
    for (i = 0; i < npar; i++) {
	eps = OS->ndeps[i]/(OS->parscale[i]);
	dpar[i] = dpar[i] + eps;
	fmingr(npar, dpar, df1, (void *)OS);
	dpar[i] = dpar[i] - 2 * eps;
	fmingr(npar, dpar, df2, (void *)OS);
	for (j = 0; j < npar; j++)
	    REAL(ans)[i * npar + j] = (OS->fnscale) * (df1[j] - df2[j])/
		(2 * eps * (OS->parscale[i]) * (OS->parscale[j]));
	dpar[i] = dpar[i] + eps;
    }
    vmaxset(vmax);
    UNPROTECT(4);
    return ans;
}
 
 
static double ** matrix(int nrh, int nch)
{
    int   i;
    double **m;
 
    m = (double **) R_alloc((nrh + 1), sizeof(double *));
    for (i = 0; i <= nrh; i++)
	m[i] = (double*) R_alloc((nch + 1), sizeof(double));
    return m;
}
 
static double ** Lmatrix(int n)
{
    int   i;
    double **m;
 
    m = (double **) R_alloc(n, sizeof(double *));
    for (i = 0; i < n; i++)
	m[i] = (double *) R_alloc((i + 1), sizeof(double));
    return m;
}
 
 
 
#define stepredn	0.2
#define acctol		0.0001
#define reltest		10.0
 
 
/*  BFGS variable-metric method, based on Pascal code
in J.C. Nash, `Compact Numerical Methods for Computers', 2nd edition,
converted by p2c then re-crafted by B.D. Ripley */
 
void
vmmin(int n0, double *b, double *Fmin, optimfn fminfn, optimgr fmingr,
      int maxit, int trace, int *mask,
      double abstol, double reltol, int nREPORT, void *ex,
      int *fncount, int *grcount, int *fail)
{
    Rboolean accpoint, enough;
    double *g, *t, *X, *c, **B;
    int   count, funcount, gradcount;
    double f, gradproj;
    int   i, j, ilast, iter = 0;
    double s, steplength;
    double D1, D2;
    int   n, *l;
 
    if (maxit <= 0) {
	*fail = 0;
	*Fmin = fminfn(n0, b, ex);
	*fncount = *grcount = 0;
	return;
    }
 
    if (nREPORT <= 0)
	error("REPORT must be > 0 (method = \"BFGS\")");
    l = (int *) R_alloc(n0, sizeof(int));
    n = 0;
    for (i = 0; i < n0; i++) if (mask[i]) l[n++] = i;
    g = vect(n0);
    t = vect(n);
    X = vect(n);
    c = vect(n);
    B = Lmatrix(n);
    f = fminfn(n0, b, ex);
    if (!R_FINITE(f))
	error("initial value in vmmin is not finite");
    if (trace) Rprintf("initial  value %f \n", f);
    *Fmin = f;
    funcount = gradcount = 1;
    fmingr(n0, b, g, ex);
    iter++;
    ilast = gradcount;
 
    do {
	if (ilast == gradcount) {
	    for (i = 0; i < n; i++) {
		for (j = 0; j < i; j++) B[i][j] = 0.0;
		B[i][i] = 1.0;
	    }
	}
	for (i = 0; i < n; i++) {
	    X[i] = b[l[i]];
	    c[i] = g[l[i]];
	}
	gradproj = 0.0;
	for (i = 0; i < n; i++) {
	    s = 0.0;
	    for (j = 0; j <= i; j++) s -= B[i][j] * g[l[j]];
	    for (j = i + 1; j < n; j++) s -= B[j][i] * g[l[j]];
	    t[i] = s;
	    gradproj += s * g[l[i]];
	}
 
	if (gradproj < 0.0) {	/* search direction is downhill */
	    steplength = 1.0;
	    accpoint = FALSE;
	    do {
		count = 0;
		for (i = 0; i < n; i++) {
		    b[l[i]] = X[i] + steplength * t[i];
		    if (reltest + X[i] == reltest + b[l[i]]) /* no change */
			count++;
		}
		if (count < n) {
		    f = fminfn(n0, b, ex);
		    funcount++;
		    accpoint = R_FINITE(f) &&
			(f <= *Fmin + gradproj * steplength * acctol);
		    if (!accpoint) {
			steplength *= stepredn;
		    }
		}
	    } while (!(count == n || accpoint));
	    enough = (f > abstol) &&
		fabs(f - *Fmin) > reltol * (fabs(*Fmin) + reltol);
	    /* stop if value if small or if relative change is low */
	    if (!enough) {
		count = n;
		*Fmin = f;
	    }
	    if (count < n) {/* making progress */
		*Fmin = f;
		fmingr(n0, b, g, ex);
		gradcount++;
		iter++;
		D1 = 0.0;
		for (i = 0; i < n; i++) {
		    t[i] = steplength * t[i];
		    c[i] = g[l[i]] - c[i];
		    D1 += t[i] * c[i];
		}
		if (D1 > 0) {
		    D2 = 0.0;
		    for (i = 0; i < n; i++) {
			s = 0.0;
			for (j = 0; j <= i; j++)
			    s += B[i][j] * c[j];
			for (j = i + 1; j < n; j++)
			    s += B[j][i] * c[j];
			X[i] = s;
			D2 += s * c[i];
		    }
		    D2 = 1.0 + D2 / D1;
		    for (i = 0; i < n; i++) {
			for (j = 0; j <= i; j++)
			    B[i][j] += (D2 * t[i] * t[j]
					- X[i] * t[j] - t[i] * X[j]) / D1;
		    }
		} else {	/* D1 < 0 */
		    ilast = gradcount;
		}
	    } else {	/* no progress */
		if (ilast < gradcount) {
		    count = 0;
		    ilast = gradcount;
		}
	    }
	} else {		/* uphill search */
	    count = 0;
	    if (ilast == gradcount) count = n;
	    else ilast = gradcount;
	    /* Resets unless has just been reset */
	}
	if (trace && (iter % nREPORT == 0))
	    Rprintf("iter%4d value %f\n", iter, f);
	if (iter >= maxit) break;
	if (gradcount - ilast > 2 * n)
	    ilast = gradcount;	/* periodic restart */
    } while (count != n || ilast != gradcount);
    if (trace) {
	Rprintf("final  value %f \n", *Fmin);
	if (iter < maxit) Rprintf("converged\n");
	else Rprintf("stopped after %i iterations\n", iter);
    }
    *fail = (iter < maxit) ? 0 : 1;
    *fncount = funcount;
    *grcount = gradcount;
}
 
 
#define big             1.0e+35   /*a very large number*/
 
 
/* Nelder-Mead */
void nmmin(int n, double *Bvec, double *X, double *Fmin, optimfn fminfn,
	   int *fail, double abstol, double intol, void *ex,
	   double alpha, double bet, double gamm, int trace,
	   int *fncount, int maxit)
{
    char action[50];
    int C;
    Rboolean calcvert, shrinkfail = FALSE;
    double convtol, f;
    int funcount=0, H, i, j, L=0;
    int n1=0;
    double oldsize;
    double **P;
    double size, step, temp, trystep;
    char tstr[6];
    double VH, VL, VR;
 
    if (maxit <= 0) {
	*Fmin = fminfn(n, Bvec, ex);
	*fncount = 0;
	*fail = 0;
	return;
    }
    if (trace)
	Rprintf("  Nelder-Mead direct search function minimizer\n");
    P = matrix(n, n+1);
    *fail = FALSE;
    f = fminfn(n, Bvec, ex);
    if (!R_FINITE(f)) {
	error("Function cannot be evaluated at initial parameters");
	*fail = TRUE;
    } else {
	if (trace) Rprintf("Function value for initial parameters = %f\n", f);
	funcount = 1;
	convtol = intol * (fabs(f) + intol);
	if (trace) Rprintf("  Scaled convergence tolerance is %g\n", convtol);
	n1 = n + 1;
	C = n + 2;
	P[n1 - 1][0] = f;
	for (i = 0; i < n; i++)
	    P[i][0] = Bvec[i];
 
	L = 1;
	size = 0.0;
 
	step = 0.0;
	for (i = 0; i < n; i++) {
	    if (0.1 * fabs(Bvec[i]) > step)
		step = 0.1 * fabs(Bvec[i]);
	}
	if (step == 0.0) step = 0.1;
	if (trace) Rprintf("Stepsize computed as %f\n", step);
	for (j = 2; j <= n1; j++) {
	    strcpy(action, "BUILD          ");
	    for (i = 0; i < n; i++)
		P[i][j - 1] = Bvec[i];
 
	    trystep = step;
	    while (P[j - 2][j - 1] == Bvec[j - 2]) {
		P[j - 2][j - 1] = Bvec[j - 2] + trystep;
		trystep *= 10;
	    }
	    size += trystep;
	}
	oldsize = size;
	calcvert = TRUE;
	shrinkfail = FALSE;
	do {
	    if (calcvert) {
		for (j = 0; j < n1; j++) {
		    if (j + 1 != L) {
			for (i = 0; i < n; i++)
			    Bvec[i] = P[i][j];
			f = fminfn(n, Bvec, ex);
			if (!R_FINITE(f)) f = big;
			funcount++;
			P[n1 - 1][j] = f;
		    }
		}
		calcvert = FALSE;
	    }
 
	    VL = P[n1 - 1][L - 1];
	    VH = VL;
	    H = L;
 
	    for (j = 1; j <= n1; j++) {
		if (j != L) {
		    f = P[n1 - 1][j - 1];
		    if (f < VL) {
			L = j;
			VL = f;
		    }
		    if (f > VH) {
			H = j;
			VH = f;
		    }
		}
	    }
 
	    if (VH > VL + convtol && VL > abstol) {
		sprintf(tstr, "%5d", funcount);
		if (trace) Rprintf("%s%s %f %f\n", action, tstr, VH, VL);
 
		for (i = 0; i < n; i++) {
		    temp = -P[i][H - 1];
		    for (j = 0; j < n1; j++)
			temp += P[i][j];
		    P[i][C - 1] = temp / n;
		}
		for (i = 0; i < n; i++)
		    Bvec[i] = (1.0 + alpha) * P[i][C - 1] - alpha * P[i][H - 1];
		f = fminfn(n, Bvec, ex);
		if (!R_FINITE(f)) f = big;
		funcount++;
		strcpy(action, "REFLECTION     ");
		VR = f;
		if (VR < VL) {
		    P[n1 - 1][C - 1] = f;
		    for (i = 0; i < n; i++) {
			f = gamm * Bvec[i] + (1 - gamm) * P[i][C - 1];
			P[i][C - 1] = Bvec[i];
			Bvec[i] = f;
		    }
		    f = fminfn(n, Bvec, ex);
		    if (!R_FINITE(f)) f = big;
		    funcount++;
		    if (f < VR) {
			for (i = 0; i < n; i++)
			    P[i][H - 1] = Bvec[i];
			P[n1 - 1][H - 1] = f;
			strcpy(action, "EXTENSION      ");
		    } else {
			for (i = 0; i < n; i++)
			    P[i][H - 1] = P[i][C - 1];
			P[n1 - 1][H - 1] = VR;
		    }
		} else {
		    strcpy(action, "HI-REDUCTION   ");
		    if (VR < VH) {
			for (i = 0; i < n; i++)
			    P[i][H - 1] = Bvec[i];
			P[n1 - 1][H - 1] = VR;
			strcpy(action, "LO-REDUCTION   ");
		    }
 
		    for (i = 0; i < n; i++)
			Bvec[i] = (1 - bet) * P[i][H - 1] + bet * P[i][C - 1];
		    f = fminfn(n, Bvec, ex);
		    if (!R_FINITE(f)) f = big;
		    funcount++;
 
		    if (f < P[n1 - 1][H - 1]) {
			for (i = 0; i < n; i++)
			    P[i][H - 1] = Bvec[i];
			P[n1 - 1][H - 1] = f;
		    } else {
			if (VR >= VH) {
			    strcpy(action, "SHRINK         ");
			    calcvert = TRUE;
			    size = 0.0;
			    for (j = 0; j < n1; j++) {
				if (j + 1 != L) {
				    for (i = 0; i < n; i++) {
					P[i][j] = bet * (P[i][j] - P[i][L - 1])
					    + P[i][L - 1];
					size += fabs(P[i][j] - P[i][L - 1]);
				    }
				}
			    }
			    if (size < oldsize) {
				shrinkfail = FALSE;
				oldsize = size;
			    } else {
				if (trace)
				    Rprintf("Polytope size measure not decreased in shrink\n");
				shrinkfail = TRUE;
			    }
			}
		    }
		}
	    }
 
	} while (!(VH <= VL + convtol || VL <= abstol ||
		   shrinkfail || funcount > maxit));
 
    }
 
    if (trace) {
	Rprintf("Exiting from Nelder Mead minimizer\n");
	Rprintf("    %d function evaluations used\n", funcount);
    }
    *Fmin = P[n1 - 1][L - 1];
    for (i = 0; i < n; i++) X[i] = P[i][L - 1];
    if (shrinkfail) *fail = 10;
    if (funcount > maxit) *fail = 1;
    *fncount = funcount;
}
 
void cgmin(int n, double *Bvec, double *X, double *Fmin,
	   optimfn fminfn, optimgr fmingr, int *fail,
	   double abstol, double intol, void *ex, int type, int trace,
	   int *fncount, int *grcount, int maxit)
{
    Rboolean accpoint;
    double *c, *g, *t;
    int count, cycle, cyclimit;
    double f;
    double G1, G2, G3, gradproj;
    int funcount=0, gradcount=0, i;
    double newstep, oldstep, setstep, steplength=1.0;
    double tol;
 
    if (maxit <= 0) {
	*Fmin = fminfn(n, Bvec, ex);
	*fncount = *grcount = 0;
	*fail = FALSE;
	return;
    }
    if (trace) {
	Rprintf("  Conjugate gradients function minimiser\n");
	switch (type) {
	case 1:	    Rprintf("Method: Fletcher Reeves\n");	break;
	case 2:	    Rprintf("Method: Polak Ribiere\n");		break;
	case 3:	    Rprintf("Method: Beale Sorenson\n");	break;
	default:
	    error("unknown type in CG method of optim");
	}
    }
    c = vect(n); g = vect(n); t = vect(n);
 
    setstep = 1.7;
    *fail = 0;
    cyclimit = n;
    tol = intol * n * sqrt(intol);
 
    if (trace) Rprintf("tolerance used in gradient test=%g\n", tol);
    f = fminfn(n, Bvec, ex);
    if (!R_FINITE(f)) {
	error("Function cannot be evaluated at initial parameters");
    } else {
	*Fmin = f;
	funcount = 1;
	gradcount = 0;
	do {
	    for (i = 0; i < n; i++) {
		t[i] = 0.0;
		c[i] = 0.0;
	    }
	    cycle = 0;
	    oldstep = 1.0;
	    count = 0;
	    do {
		cycle++;
		if (trace) {
		    Rprintf("%d %d %f\n", gradcount, funcount, *Fmin);
		    Rprintf("parameters ");
		    for (i = 1; i <= n; i++) {
			Rprintf("%10.5f ", Bvec[i - 1]);
			if (i / 7 * 7 == i && i < n)
			    Rprintf("\n");
		    }
		    Rprintf("\n");
		}
		gradcount++;
		if (gradcount > maxit) {
		    *fncount = funcount;
		    *grcount = gradcount;
		    *fail = 1;
		    return;
		}
		fmingr(n, Bvec, g, ex);
		G1 = 0.0;
		G2 = 0.0;
		for (i = 0; i < n; i++) {
		    X[i] = Bvec[i];
		    switch (type) {
 
		    case 1: /* Fletcher-Reeves */
			G1 += g[i] * g[i];
			G2 += c[i] * c[i];
			break;
 
		    case 2: /* Polak-Ribiere */
			G1 += g[i] * (g[i] - c[i]);
			G2 += c[i] * c[i];
			break;
 
		    case 3: /* Beale-Sorenson */
			G1 += g[i] * (g[i] - c[i]);
			G2 += t[i] * (g[i] - c[i]);
			break;
 
		    default:
			error("unknown type in CG method of optim");
		    }
		    c[i] = g[i];
		}
		if (G1 > tol) {
		    if (G2 > 0.0)
			G3 = G1 / G2;
		    else
			G3 = 1.0;
		    gradproj = 0.0;
		    for (i = 0; i < n; i++) {
			t[i] = t[i] * G3 - g[i];
			gradproj += t[i] * g[i];
		    }
		    steplength = oldstep;
 
		    accpoint = FALSE;
		    do {
			count = 0;
			for (i = 0; i < n; i++) {
			    Bvec[i] = X[i] + steplength * t[i];
			    if (reltest + X[i] == reltest + Bvec[i])
				count++;
			}
			if (count < n) {
			    f = fminfn(n, Bvec, ex);
			    funcount++;
			    accpoint = (R_FINITE(f) &&
					f <= *Fmin + gradproj * steplength * acctol);
 
			    if (!accpoint) {
				steplength *= stepredn;
				if (trace) Rprintf("*");
			    }
			}
		    } while (!(count == n || accpoint));
		    if (count < n) {
			newstep = 2 * (f - *Fmin - gradproj * steplength);
			if (newstep > 0) {
			    newstep = -(gradproj * steplength * steplength / newstep);
			    for (i = 0; i < n; i++)
				Bvec[i] = X[i] + newstep * t[i];
			    *Fmin = f;
			    f = fminfn(n, Bvec, ex);
			    funcount++;
			    if (f < *Fmin) {
				*Fmin = f;
				if (trace) Rprintf(" i< ");
			    } else {
				if (trace) Rprintf(" i> ");
				for (i = 0; i < n; i++)
				    Bvec[i] = X[i] + steplength * t[i];
			    }
			}
		    }
		}
		oldstep = setstep * steplength;
		if (oldstep > 1.0)
		    oldstep = 1.0;
	    } while ((count != n) && (G1 > tol) && (cycle != cyclimit));
 
	} while ((cycle != 1) ||
		 ((count != n) && (G1 > tol) && *Fmin > abstol));
 
    }
    if (trace) {
	Rprintf("Exiting from conjugate gradients minimizer\n");
	Rprintf("    %d function evaluations used\n", funcount);
	Rprintf("    %d gradient evaluations used\n", gradcount);
    }
    *fncount = funcount;
    *grcount = gradcount;
}
 
void lbfgsb(int n, int m, double *x, double *l, double *u, int *nbd,
	    double *Fmin, optimfn fminfn, optimgr fmingr, int *fail,
	    void *ex, double factr, double pgtol,
	    int *fncount, int *grcount, int maxit, char *msg,
	    int trace, int nREPORT)
{
    char task[60];
    double f, *g, dsave[29], *wa;
    int tr = -1, iter = 0, *iwa, isave[44], lsave[4];
 
    if (nREPORT <= 0)
	error("REPORT must be > 0 (method = \"L-BFGS-B\")");
    switch(trace) {
    case 2: tr = 0; break;
    case 3: tr = nREPORT; break;
    case 4: tr = 99; break;
    case 5: tr = 100; break;
    case 6: tr = 101; break;
    default: tr = -1; break;
    }
 
    *fail = 0;
    g = vect(n);
    wa = vect(2*m*n+4*n+11*m*m+8*m);
    iwa = (int *) R_alloc(3*n, sizeof(int));
    strcpy(task, "START");
    while(1) {
	/* Main workhorse setulb() from ../appl/lbfgsb.c : */
	setulb(n, m, x, l, u, nbd, &f, g, factr, &pgtol, wa, iwa, task,
	       tr, lsave, isave, dsave);
/*	Rprintf("in lbfgsb - %s\n", task);*/
	if (strncmp(task, "FG", 2) == 0) {
	    f = fminfn(n, x, ex);
	    if (!R_FINITE(f))
		error("L-BFGS-B needs finite values of fn");
	    fmingr(n, x, g, ex);
	} else if (strncmp(task, "NEW_X", 5) == 0) {
	    if(trace == 1 && (iter % nREPORT == 0)) {
		Rprintf("iter %4d value %f\n", iter, f);
	    }
	    if (++iter > maxit) {
		*fail = 1;
		break;
	    }
	} else if (strncmp(task, "WARN", 4) == 0) {
	    *fail = 51;
	    break;
	} else if (strncmp(task, "CONV", 4) == 0) {
	    break;
	} else if (strncmp(task, "ERROR", 5) == 0) {
	    *fail = 52;
	    break;
	} else { /* some other condition that is not supposed to happen */
	    *fail = 52;
	    break;
	}
    }
    *Fmin = f;
    *fncount = *grcount = isave[33];
    if (trace) {
	Rprintf("final  value %f \n", *Fmin);
	if (iter < maxit && *fail == 0) Rprintf("converged\n");
	else Rprintf("stopped after %i iterations\n", iter);
    }
    strcpy(msg, task);
}
 
 
#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)
 
/* Given a starting point pb[0..n-1], simulated annealing minimization
   is performed on the function fminfn. The starting temperature
   is input as ti. To make sann work silently set trace to zero.
   sann makes in total maxit function evaluations, tmax
   evaluations at each temperature. Returned quantities are pb
   (the location of the minimum), and yb (the minimum value of
   the function func).  Author: Adrian Trapletti
*/
{
    long j;
    int k, its, itdoc;
    double t, y, dy, ytry, scale;
    double *p, *dp, *ptry;
 
    p = vect (n); dp = vect (n); ptry = vect (n);
    GetRNGstate();
    *yb = fminfn (n, pb, ex);  /* init best system state pb, *yb */
    if (!R_FINITE(*yb)) *yb = big;
    for (j = 0; j < n; j++) p[j] = pb[j];
    y = *yb;  /* init system state p, y */
    if (trace)
    {
	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];
		y = ytry;  /* update system state p, y */
		if (y <= *yb)  /* if system state is best, then update best system state pb, *yb */
		{
		    for (j = 0; j < n; j++) pb[j] = p[j];
		    *yb = y;
		}
	    }
	    its++; k++;
	}
	if ((trace) && ((itdoc % STEPS) == 0))
	    Rprintf("iter %8d value %f\n", its - 1, *yb);
	itdoc++;
    }
    if (trace)
    {
	Rprintf ("final         value %f\n", *yb);
	Rprintf ("sann stopped after %d iterations\n", its - 1);
    }
    PutRNGstate();
}
 
#undef E1
#undef STEPS

Rev 19233	Rev 21689
1	`/*`	1	`/*`
2	`* R : A Computer Language for Statistical Data Analysis`	2	`* R : A Computer Language for Statistical Data Analysis`
3	`* Copyright (C) 1999-2002 the R Development Core Team`	3	`* Copyright (C) 1999-2002 the R Development Core Team`
4	`*`	4	`*`
5	`* This program is free software; you can redistribute it and/or modify`	5	`* This program is free software; you can redistribute it and/or modify`
6	`* it under the terms of the GNU General Public License as published by`	6	`* it under the terms of the GNU General Public License as published by`
7	`* the Free Software Foundation; either version 2 of the License, or`	7	`* the Free Software Foundation; either version 2 of the License, or`
8	`* (at your option) any later version.`	8	`* (at your option) any later version.`
9	`*`	9	`*`
10	`* This program is distributed in the hope that it will be useful,`	10	`* This program is distributed in the hope that it will be useful,`
11	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`	11	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`
12	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`	12	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`
13	`* GNU General Public License for more details.`	13	`* GNU General Public License for more details.`
14	`*`	14	`*`
15	`* You should have received a copy of the GNU General Public License`	15	`* You should have received a copy of the GNU General Public License`
16	`* along with this program; if not, write to the Free Software`	16	`* along with this program; if not, write to the Free Software`
17	`* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA`	17	`* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA`
18	`*/`	18	`*/`
19		19
20	`#include <Defn.h>`	20	`#include <Defn.h>`
21	`#include <R_ext/Random.h> /* for the random number generation in`	21	`#include <R_ext/Random.h> /* for the random number generation in`
22	`samin() */`	22	`samin() */`
23	`#include <R_ext/Applic.h> /* setulb() */`	23	`#include <R_ext/Applic.h> /* setulb() */`
24		24
25	`static SEXP getListElement(SEXP list, char *str)`	25	`static SEXP getListElement(SEXP list, char *str)`
26	`{`	26	`{`
27	`SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);`	27	`SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);`
28	`int i;`	28	`int i;`
29		29
30	`for (i = 0; i < length(list); i++)`	30	`for (i = 0; i < length(list); i++)`
31	`if (strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {`	31	`if (strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {`
32	`elmt = VECTOR_ELT(list, i);`	32	`elmt = VECTOR_ELT(list, i);`
33	`break;`	33	`break;`
34	`}`	34	`}`
35	`return elmt;`	35	`return elmt;`
36	`}`	36	`}`
37		37
38	`static double * vect(int n)`	38	`static double * vect(int n)`
39	`{`	39	`{`
40	`return (double *)R_alloc(n, sizeof(double));`	40	`return (double *)R_alloc(n, sizeof(double));`
41	`}`	41	`}`
42		42
43	`typedef struct opt_struct`	43	`typedef struct opt_struct`
44	`{`	44	`{`
45	`SEXP R_fcall; /* function */`	45	`SEXP R_fcall; /* function */`
46	`SEXP R_gcall; /* gradient */`	46	`SEXP R_gcall; /* gradient */`
47	`SEXP R_env; /* where to evaluate the calls */`	47	`SEXP R_env; /* where to evaluate the calls */`
48	`double* ndeps; /* tolerances for numerical derivatives */`	48	`double* ndeps; /* tolerances for numerical derivatives */`
49	`double fnscale; /* scaling for objective */`	49	`double fnscale; /* scaling for objective */`
50	`double* parscale;/* scaling for parameters */`	50	`double* parscale;/* scaling for parameters */`
51	`int usebounds;`	51	`int usebounds;`
52	`double* lower, *upper;`	52	`double* lower, *upper;`
53	`} opt_struct, *OptStruct;`	53	`} opt_struct, *OptStruct;`
54		54
55		55
56		56
57	`static double fminfn(int n, double p, void ex)`	57	`static double fminfn(int n, double p, void ex)`
58	`{`	58	`{`
59	`SEXP s, x;`	59	`SEXP s, x;`
60	`int i;`	60	`int i;`
61	`double val;`	61	`double val;`
62	`OptStruct OS = (OptStruct) ex;`	62	`OptStruct OS = (OptStruct) ex;`
63	`PROTECT_INDEX ipx;`	63	`PROTECT_INDEX ipx;`
64		64
65	`PROTECT(x = allocVector(REALSXP, n));`	65	`PROTECT(x = allocVector(REALSXP, n));`
66	`for (i = 0; i < n; i++) {`	66	`for (i = 0; i < n; i++) {`
67	`if (!R_FINITE(p[i])) error("non-finite value supplied by optim");`	67	`if (!R_FINITE(p[i])) error("non-finite value supplied by optim");`
68	`REAL(x)[i] = p[i] * (OS->parscale[i]);`	68	`REAL(x)[i] = p[i] * (OS->parscale[i]);`
69	`}`	69	`}`
70	`SETCADR(OS->R_fcall, x);`	70	`SETCADR(OS->R_fcall, x);`
71	`PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);`	71	`PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);`
72	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`	72	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
73	`val = REAL(s)[0]/(OS->fnscale);`	73	`val = REAL(s)[0]/(OS->fnscale);`
74	`UNPROTECT(2);`	74	`UNPROTECT(2);`
75	`return val;`	75	`return val;`
76	`}`	76	`}`
77		77
78	`static void fmingr(int n, double p, double df, void *ex)`	78	`static void fmingr(int n, double p, double df, void *ex)`
79	`{`	79	`{`
80	`SEXP s, x;`	80	`SEXP s, x;`
81	`int i;`	81	`int i;`
82	`double val1, val2, eps, epsused, tmp;`	82	`double val1, val2, eps, epsused, tmp;`
83	`OptStruct OS = (OptStruct) ex;`	83	`OptStruct OS = (OptStruct) ex;`
84	`PROTECT_INDEX ipx;`	84	`PROTECT_INDEX ipx;`
85		85
86	`if (!isNull(OS->R_gcall)) { /* analytical derivatives */`	86	`if (!isNull(OS->R_gcall)) { /* analytical derivatives */`
87	`PROTECT(x = allocVector(REALSXP, n));`	87	`PROTECT(x = allocVector(REALSXP, n));`
88	`for (i = 0; i < n; i++) {`	88	`for (i = 0; i < n; i++) {`
89	`if (!R_FINITE(p[i])) error("non-finite value supplied by optim");`	89	`if (!R_FINITE(p[i])) error("non-finite value supplied by optim");`
90	`REAL(x)[i] = p[i] * (OS->parscale[i]);`	90	`REAL(x)[i] = p[i] * (OS->parscale[i]);`
91	`}`	91	`}`
92	`SETCADR(OS->R_gcall, x);`	92	`SETCADR(OS->R_gcall, x);`
93	`PROTECT_WITH_INDEX(s = eval(OS->R_gcall, OS->R_env), &ipx);`	93	`PROTECT_WITH_INDEX(s = eval(OS->R_gcall, OS->R_env), &ipx);`
94	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`	94	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
95	`if(LENGTH(s) != n)`	95	`if(LENGTH(s) != n)`
96	`error("gradient in optim evaluated to length %d not %d",`	96	`error("gradient in optim evaluated to length %d not %d",`
97	`LENGTH(s), n);`	97	`LENGTH(s), n);`
98	`for (i = 0; i < n; i++)`	98	`for (i = 0; i < n; i++)`
99	`df[i] = REAL(s)[i] * (OS->parscale[i])/(OS->fnscale);`	99	`df[i] = REAL(s)[i] * (OS->parscale[i])/(OS->fnscale);`
100	`UNPROTECT(2);`	100	`UNPROTECT(2);`
101	`} else { /* numerical derivatives */`	101	`} else { /* numerical derivatives */`
102	`PROTECT(x = allocVector(REALSXP, n));`	102	`PROTECT(x = allocVector(REALSXP, n));`
103	`for (i = 0; i < n; i++) REAL(x)[i] = p[i] * (OS->parscale[i]);`	103	`for (i = 0; i < n; i++) REAL(x)[i] = p[i] * (OS->parscale[i]);`
104	`SETCADR(OS->R_fcall, x);`	104	`SETCADR(OS->R_fcall, x);`
105	`if(OS->usebounds == 0) {`	105	`if(OS->usebounds == 0) {`
106	`for (i = 0; i < n; i++) {`	106	`for (i = 0; i < n; i++) {`
107	`eps = OS->ndeps[i];`	107	`eps = OS->ndeps[i];`
108	`REAL(x)[i] = (p[i] + eps) * (OS->parscale[i]);`	108	`REAL(x)[i] = (p[i] + eps) * (OS->parscale[i]);`
109	`SETCADR(OS->R_fcall, x);`	109	`SETCADR(OS->R_fcall, x);`
110	`PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);`	110	`PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);`
111	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`	111	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
112	`val1 = REAL(s)[0]/(OS->fnscale);`	112	`val1 = REAL(s)[0]/(OS->fnscale);`
113	`REAL(x)[i] = (p[i] - eps) * (OS->parscale[i]);`	113	`REAL(x)[i] = (p[i] - eps) * (OS->parscale[i]);`
114	`SETCADR(OS->R_fcall, x);`	114	`SETCADR(OS->R_fcall, x);`
115	`REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);`	115	`REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);`
116	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`	116	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
117	`val2 = REAL(s)[0]/(OS->fnscale);`	117	`val2 = REAL(s)[0]/(OS->fnscale);`
118	`df[i] = (val1 - val2)/(2 * eps);`	118	`df[i] = (val1 - val2)/(2 * eps);`
119	`#define DO_df_x \`	119	`#define DO_df_x \`
120	`if(!R_FINITE(df[i])) \`	120	`if(!R_FINITE(df[i])) \`
121	`error("non-finite finite-difference value [%d]", i);\`	121	`error("non-finite finite-difference value [%d]", i);\`
122	`REAL(x)[i] = p[i] * (OS->parscale[i])`	122	`REAL(x)[i] = p[i] * (OS->parscale[i])`
123		123
124	`DO_df_x;`	124	`DO_df_x;`
125	`UNPROTECT(1);`	125	`UNPROTECT(1);`
126	`}`	126	`}`
127	`} else { /* usebounds */`	127	`} else { /* usebounds */`
128	`for (i = 0; i < n; i++) {`	128	`for (i = 0; i < n; i++) {`
129	`epsused = eps = OS->ndeps[i];`	129	`epsused = eps = OS->ndeps[i];`
130	`tmp = p[i] + eps;`	130	`tmp = p[i] + eps;`
131	`if (tmp > OS->upper[i]) {`	131	`if (tmp > OS->upper[i]) {`
132	`tmp = OS->upper[i];`	132	`tmp = OS->upper[i];`
133	`epsused = tmp - p[i] ;`	133	`epsused = tmp - p[i] ;`
134	`}`	134	`}`
135	`REAL(x)[i] = tmp * (OS->parscale[i]);`	135	`REAL(x)[i] = tmp * (OS->parscale[i]);`
136	`SETCADR(OS->R_fcall, x);`	136	`SETCADR(OS->R_fcall, x);`
137	`PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);`	137	`PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);`
138	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`	138	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
139	`val1 = REAL(s)[0]/(OS->fnscale);`	139	`val1 = REAL(s)[0]/(OS->fnscale);`
140	`tmp = p[i] - eps;`	140	`tmp = p[i] - eps;`
141	`if (tmp < OS->lower[i]) {`	141	`if (tmp < OS->lower[i]) {`
142	`tmp = OS->lower[i];`	142	`tmp = OS->lower[i];`
143	`eps = p[i] - tmp;`	143	`eps = p[i] - tmp;`
144	`}`	144	`}`
145	`REAL(x)[i] = tmp * (OS->parscale[i]);`	145	`REAL(x)[i] = tmp * (OS->parscale[i]);`
146	`SETCADR(OS->R_fcall, x);`	146	`SETCADR(OS->R_fcall, x);`
147	`REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);`	147	`REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);`
148	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`	148	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
149	`val2 = REAL(s)[0]/(OS->fnscale);`	149	`val2 = REAL(s)[0]/(OS->fnscale);`
150	`df[i] = (val1 - val2)/(epsused + eps);`	150	`df[i] = (val1 - val2)/(epsused + eps);`
151		151
152	`DO_df_x;`	152	`DO_df_x;`
153	`UNPROTECT(1);`	153	`UNPROTECT(1);`
154	`}`	154	`}`
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;`
165	`int i, npar=0, *mask, trace, maxit, fncount, grcount, nREPORT, tmax;`	195	`int i, npar=0, *mask, trace, maxit, fncount, grcount, nREPORT, tmax;`
166	`int ifail = 0;`	196	`int ifail = 0;`
167	`double dpar, opar, val, abstol, reltol, temp;`	197	`double dpar, opar, val, abstol, reltol, temp;`
168	`char *tn;`	198	`char *tn;`
169	`OptStruct OS;`	199	`OptStruct OS;`
170	`char *vmax;`	200	`char *vmax;`
171		201
172	`checkArity(op, args);`	202	`checkArity(op, args);`
173	`vmax = vmaxget();`	203	`vmax = vmaxget();`
174	`OS = (OptStruct) R_alloc(1, sizeof(opt_struct));`	204	`OS = (OptStruct) R_alloc(1, sizeof(opt_struct));`
175	`OS->usebounds = 0;`	205	`OS->usebounds = 0;`
176	`OS->R_env = rho;`	206	`OS->R_env = rho;`
177	`par = CAR(args);`	207	`par = CAR(args);`
178	`args = CDR(args); fn = CAR(args);`	208	`args = CDR(args); fn = CAR(args);`
179	`if (!isFunction(fn)) errorcall(call, "fn is not a function");`	209	`if (!isFunction(fn)) errorcall(call, "fn is not a function");`
180	`args = CDR(args); gr = CAR(args);`	210	`args = CDR(args); gr = CAR(args);`
181	`args = CDR(args); method = CAR(args);`	211	`args = CDR(args); method = CAR(args);`
182	`if (!isString(method)\|\| LENGTH(method) != 1)`	212	`if (!isString(method)\|\| LENGTH(method) != 1)`
183	`errorcall(call, "invalid method argument");`	213	`errorcall(call, "invalid method argument");`
184	`tn = CHAR(STRING_ELT(method, 0));`	214	`tn = CHAR(STRING_ELT(method, 0));`
185	`args = CDR(args); options = CAR(args);`	215	`args = CDR(args); options = CAR(args);`
186	`PROTECT(OS->R_fcall = lang2(fn, R_NilValue));`	216	`PROTECT(OS->R_fcall = lang2(fn, R_NilValue));`
187	`PROTECT(par = coerceVector(duplicate(par), REALSXP));`	217	`PROTECT(par = coerceVector(duplicate(par), REALSXP));`
188	`npar = LENGTH(par);`	218	`npar = LENGTH(par);`
189	`dpar = vect(npar);`	219	`dpar = vect(npar);`
190	`opar = vect(npar);`	220	`opar = vect(npar);`
191	`trace = asInteger(getListElement(options, "trace"));`	221	`trace = asInteger(getListElement(options, "trace"));`
192	`OS->fnscale = asReal(getListElement(options, "fnscale"));`	222	`OS->fnscale = asReal(getListElement(options, "fnscale"));`
193	`tmp = getListElement(options, "parscale");`	223	`tmp = getListElement(options, "parscale");`
194	`if (LENGTH(tmp) != npar)`	224	`if (LENGTH(tmp) != npar)`
195	`errorcall(call, "parscale is of the wrong length");`	225	`errorcall(call, "parscale is of the wrong length");`
196	`PROTECT(tmp = coerceVector(tmp, REALSXP));`	226	`PROTECT(tmp = coerceVector(tmp, REALSXP));`
197	`OS->parscale = vect(npar);`	227	`OS->parscale = vect(npar);`
198	`for (i = 0; i < npar; i++) OS->parscale[i] = REAL(tmp)[i];`	228	`for (i = 0; i < npar; i++) OS->parscale[i] = REAL(tmp)[i];`
199	`UNPROTECT(1);`	229	`UNPROTECT(1);`
200	`for (i = 0; i < npar; i++)`	230	`for (i = 0; i < npar; i++)`
201	`dpar[i] = REAL(par)[i] / (OS->parscale[i]);`	231	`dpar[i] = REAL(par)[i] / (OS->parscale[i]);`
202	`PROTECT(res = allocVector(VECSXP, 5));`	232	`PROTECT(res = allocVector(VECSXP, 5));`
203	`PROTECT(value = allocVector(REALSXP, 1));`	233	`PROTECT(value = allocVector(REALSXP, 1));`
204	`PROTECT(counts = allocVector(INTSXP, 2));`	234	`PROTECT(counts = allocVector(INTSXP, 2));`
205	`PROTECT(conv = allocVector(INTSXP, 1));`	235	`PROTECT(conv = allocVector(INTSXP, 1));`
206	`abstol = asReal(getListElement(options, "abstol"));`	236	`abstol = asReal(getListElement(options, "abstol"));`
207	`reltol = asReal(getListElement(options, "reltol"));`	237	`reltol = asReal(getListElement(options, "reltol"));`
208	`maxit = asInteger(getListElement(options, "maxit"));`	238	`maxit = asInteger(getListElement(options, "maxit"));`
209	`if (maxit == NA_INTEGER) error("maxit is not an integer");`	239	`if (maxit == NA_INTEGER) error("maxit is not an integer");`
210		240
211	`if (strcmp(tn, "Nelder-Mead") == 0) {`	241	`if (strcmp(tn, "Nelder-Mead") == 0) {`
212	`double alpha, beta, gamm;`	242	`double alpha, beta, gamm;`
213		243
214	`alpha = asReal(getListElement(options, "alpha"));`	244	`alpha = asReal(getListElement(options, "alpha"));`
215	`beta = asReal(getListElement(options, "beta"));`	245	`beta = asReal(getListElement(options, "beta"));`
216	`gamm = asReal(getListElement(options, "gamma"));`	246	`gamm = asReal(getListElement(options, "gamma"));`
217	`nmmin(npar, dpar, opar, &val, fminfn, &ifail, abstol, reltol,`	247	`nmmin(npar, dpar, opar, &val, fminfn, &ifail, abstol, reltol,`
218	`(void *)OS, alpha, beta, gamm, trace, &fncount, maxit);`	248	`(void *)OS, alpha, beta, gamm, trace, &fncount, maxit);`
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"));`
238	`if (!isNull(gr)) {`	275	`if (!isNull(gr)) {`
239	`if (!isFunction(gr)) error("gr is not a function");`	276	`if (!isFunction(gr)) error("gr is not a function");`
240	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`	277	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`
241	`} else {`	278	`} else {`
242	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`	279	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`
243	`ndeps = getListElement(options, "ndeps");`	280	`ndeps = getListElement(options, "ndeps");`
244	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`	281	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`
245	`OS->ndeps = vect(npar);`	282	`OS->ndeps = vect(npar);`
246	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`	283	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`
247	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`	284	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`
248	`UNPROTECT(1);`	285	`UNPROTECT(1);`
249	`}`	286	`}`
250	`mask = (int *) R_alloc(npar, sizeof(int));`	287	`mask = (int *) R_alloc(npar, sizeof(int));`
251	`for (i = 0; i < npar; i++) mask[i] = 1;`	288	`for (i = 0; i < npar; i++) mask[i] = 1;`
252	`vmmin(npar, dpar, &val, fminfn, fmingr, maxit, trace, mask, abstol,`	289	`vmmin(npar, dpar, &val, fminfn, fmingr, maxit, trace, mask, abstol,`
253	`reltol, nREPORT, (void *)OS, &fncount, &grcount, &ifail);`	290	`reltol, nREPORT, (void *)OS, &fncount, &grcount, &ifail);`
254	`for (i = 0; i < npar; i++)`	291	`for (i = 0; i < npar; i++)`
255	`REAL(par)[i] = dpar[i] * (OS->parscale[i]);`	292	`REAL(par)[i] = dpar[i] * (OS->parscale[i]);`
256	`UNPROTECT(1); /* OS->R_gcall */`	293	`UNPROTECT(1); /* OS->R_gcall */`
257	`} else if (strcmp(tn, "CG") == 0) {`	294	`} else if (strcmp(tn, "CG") == 0) {`
258	`int type;`	295	`int type;`
259	`SEXP ndeps;`	296	`SEXP ndeps;`
260		297
261	`type = asInteger(getListElement(options, "type"));`	298	`type = asInteger(getListElement(options, "type"));`
262	`if (!isNull(gr)) {`	299	`if (!isNull(gr)) {`
263	`if (!isFunction(gr)) error("gr is not a function");`	300	`if (!isFunction(gr)) error("gr is not a function");`
264	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`	301	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`
265	`} else {`	302	`} else {`
266	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`	303	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`
267	`ndeps = getListElement(options, "ndeps");`	304	`ndeps = getListElement(options, "ndeps");`
268	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`	305	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`
269	`OS->ndeps = vect(npar);`	306	`OS->ndeps = vect(npar);`
270	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`	307	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`
271	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`	308	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`
272	`UNPROTECT(1);`	309	`UNPROTECT(1);`
273	`}`	310	`}`
274	`cgmin(npar, dpar, opar, &val, fminfn, fmingr, &ifail, abstol,`	311	`cgmin(npar, dpar, opar, &val, fminfn, fmingr, &ifail, abstol,`
275	`reltol, (void *)OS, type, trace, &fncount, &grcount, maxit);`	312	`reltol, (void *)OS, type, trace, &fncount, &grcount, maxit);`
276	`for (i = 0; i < npar; i++)`	313	`for (i = 0; i < npar; i++)`
277	`REAL(par)[i] = opar[i] * (OS->parscale[i]);`	314	`REAL(par)[i] = opar[i] * (OS->parscale[i]);`
278	`UNPROTECT(1); /* OS->R_gcall */`	315	`UNPROTECT(1); /* OS->R_gcall */`
279		316
280	`} else if (strcmp(tn, "L-BFGS-B") == 0) {`	317	`} else if (strcmp(tn, "L-BFGS-B") == 0) {`
281	`SEXP ndeps, smsg;`	318	`SEXP ndeps, smsg;`
282	`double lower = vect(npar), upper = vect(npar);`	319	`double lower = vect(npar), upper = vect(npar);`
283	`int lmm, nbd = (int ) R_alloc(npar, sizeof(int));`	320	`int lmm, nbd = (int ) R_alloc(npar, sizeof(int));`
284	`double factr, pgtol;`	321	`double factr, pgtol;`
285	`char msg[60];`	322	`char msg[60];`
286		323
287	`nREPORT = asInteger(getListElement(options, "REPORT"));`	324	`nREPORT = asInteger(getListElement(options, "REPORT"));`
288	`factr = asReal(getListElement(options, "factr"));`	325	`factr = asReal(getListElement(options, "factr"));`
289	`pgtol = asReal(getListElement(options, "pgtol"));`	326	`pgtol = asReal(getListElement(options, "pgtol"));`
290	`lmm = asInteger(getListElement(options, "lmm"));`	327	`lmm = asInteger(getListElement(options, "lmm"));`
291	`if (!isNull(gr)) {`	328	`if (!isNull(gr)) {`
292	`if (!isFunction(gr)) error("gr is not a function");`	329	`if (!isFunction(gr)) error("gr is not a function");`
293	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`	330	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`
294	`} else {`	331	`} else {`
295	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`	332	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`
296	`ndeps = getListElement(options, "ndeps");`	333	`ndeps = getListElement(options, "ndeps");`
297	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`	334	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`
298	`OS->ndeps = vect(npar);`	335	`OS->ndeps = vect(npar);`
299	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`	336	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`
300	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`	337	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`
301	`UNPROTECT(1);`	338	`UNPROTECT(1);`
302	`}`	339	`}`
303	`args = CDR(args); slower = CAR(args); /* coerce in calling code */`	340	`args = CDR(args); slower = CAR(args); /* coerce in calling code */`
304	`args = CDR(args); supper = CAR(args);`	341	`args = CDR(args); supper = CAR(args);`
305	`for (i = 0; i < npar; i++) {`	342	`for (i = 0; i < npar; i++) {`
306	`lower[i] = REAL(slower)[i] / (OS->parscale[i]);`	343	`lower[i] = REAL(slower)[i] / (OS->parscale[i]);`
307	`upper[i] = REAL(supper)[i] / (OS->parscale[i]);`	344	`upper[i] = REAL(supper)[i] / (OS->parscale[i]);`
308	`if (!R_FINITE(lower[i])) {`	345	`if (!R_FINITE(lower[i])) {`
309	`if (!R_FINITE(upper[i])) nbd[i] = 0; else nbd[i] = 3;`	346	`if (!R_FINITE(upper[i])) nbd[i] = 0; else nbd[i] = 3;`
310	`} else {`	347	`} else {`
311	`if (!R_FINITE(upper[i])) nbd[i] = 1; else nbd[i] = 2;`	348	`if (!R_FINITE(upper[i])) nbd[i] = 1; else nbd[i] = 2;`
312	`}`	349	`}`
313	`}`	350	`}`
314	`OS->usebounds = 1;`	351	`OS->usebounds = 1;`
315	`OS->lower = lower;`	352	`OS->lower = lower;`
316	`OS->upper = upper;`	353	`OS->upper = upper;`
317	`lbfgsb(npar, lmm, dpar, lower, upper, nbd, &val, fminfn, fmingr,`	354	`lbfgsb(npar, lmm, dpar, lower, upper, nbd, &val, fminfn, fmingr,`
318	`&ifail, (void *)OS, factr, pgtol, &fncount, &grcount,`	355	`&ifail, (void *)OS, factr, pgtol, &fncount, &grcount,`
319	`maxit, msg, trace, nREPORT);`	356	`maxit, msg, trace, nREPORT);`
320	`for (i = 0; i < npar; i++)`	357	`for (i = 0; i < npar; i++)`
321	`REAL(par)[i] = dpar[i] * (OS->parscale[i]);`	358	`REAL(par)[i] = dpar[i] * (OS->parscale[i]);`
322	`UNPROTECT(1); /* OS->R_gcall */`	359	`UNPROTECT(1); /* OS->R_gcall */`
323	`PROTECT(smsg = allocVector(STRSXP, 1));`	360	`PROTECT(smsg = allocVector(STRSXP, 1));`
324	`SET_STRING_ELT(smsg, 0, mkChar(msg));`	361	`SET_STRING_ELT(smsg, 0, mkChar(msg));`
325	`SET_VECTOR_ELT(res, 4, smsg);`	362	`SET_VECTOR_ELT(res, 4, smsg);`
326	`UNPROTECT(1);`	363	`UNPROTECT(1);`
327	`} else`	364	`} else`
328	`errorcall(call, "unknown method");`	365	`errorcall(call, "unknown method");`
329		366
330	`REAL(value)[0] = val * (OS->fnscale);`	367	`REAL(value)[0] = val * (OS->fnscale);`
331	`SET_VECTOR_ELT(res, 0, par); SET_VECTOR_ELT(res, 1, value);`	368	`SET_VECTOR_ELT(res, 0, par); SET_VECTOR_ELT(res, 1, value);`
332	`INTEGER(counts)[0] = fncount; INTEGER(counts)[1] = grcount;`	369	`INTEGER(counts)[0] = fncount; INTEGER(counts)[1] = grcount;`
333	`SET_VECTOR_ELT(res, 2, counts);`	370	`SET_VECTOR_ELT(res, 2, counts);`
334	`INTEGER(conv)[0] = ifail;`	371	`INTEGER(conv)[0] = ifail;`
335	`SET_VECTOR_ELT(res, 3, conv);`	372	`SET_VECTOR_ELT(res, 3, conv);`
336	`vmaxset(vmax);`	373	`vmaxset(vmax);`
337	`UNPROTECT(6);`	374	`UNPROTECT(6);`
338	`return res;`	375	`return res;`
339	`}`	376	`}`
340		377
341	`/* par fn gr options */`	378	`/* par fn gr options */`
342	`SEXP do_optimhess(SEXP call, SEXP op, SEXP args, SEXP rho)`	379	`SEXP do_optimhess(SEXP call, SEXP op, SEXP args, SEXP rho)`
343	`{`	380	`{`
344	`SEXP par, fn, gr, options, tmp, ndeps, ans;`	381	`SEXP par, fn, gr, options, tmp, ndeps, ans;`
345	`OptStruct OS;`	382	`OptStruct OS;`
346	`int npar, i , j;`	383	`int npar, i , j;`
347	`double dpar, df1, *df2, eps;`	384	`double dpar, df1, *df2, eps;`
348	`char *vmax;`	385	`char *vmax;`
349		386
350	`checkArity(op, args);`	387	`checkArity(op, args);`
351	`vmax = vmaxget();`	388	`vmax = vmaxget();`
352	`OS = (OptStruct) R_alloc(1, sizeof(opt_struct));`	389	`OS = (OptStruct) R_alloc(1, sizeof(opt_struct));`
353	`OS->usebounds = 0;`	390	`OS->usebounds = 0;`
354	`OS->R_env = rho;`	391	`OS->R_env = rho;`
355	`par = CAR(args);`	392	`par = CAR(args);`
356	`npar = LENGTH(par);`	393	`npar = LENGTH(par);`
357	`args = CDR(args); fn = CAR(args);`	394	`args = CDR(args); fn = CAR(args);`
358	`if (!isFunction(fn)) errorcall(call, "fn is not a function");`	395	`if (!isFunction(fn)) errorcall(call, "fn is not a function");`
359	`args = CDR(args); gr = CAR(args);`	396	`args = CDR(args); gr = CAR(args);`
360	`args = CDR(args); options = CAR(args);`	397	`args = CDR(args); options = CAR(args);`
361	`OS->fnscale = asReal(getListElement(options, "fnscale"));`	398	`OS->fnscale = asReal(getListElement(options, "fnscale"));`
362	`tmp = getListElement(options, "parscale");`	399	`tmp = getListElement(options, "parscale");`
363	`if (LENGTH(tmp) != npar)`	400	`if (LENGTH(tmp) != npar)`
364	`errorcall(call, "parscale is of the wrong length");`	401	`errorcall(call, "parscale is of the wrong length");`
365	`PROTECT(tmp = coerceVector(tmp, REALSXP));`	402	`PROTECT(tmp = coerceVector(tmp, REALSXP));`
366	`OS->parscale = vect(npar);`	403	`OS->parscale = vect(npar);`
367	`for (i = 0; i < npar; i++) OS->parscale[i] = REAL(tmp)[i];`	404	`for (i = 0; i < npar; i++) OS->parscale[i] = REAL(tmp)[i];`
368	`UNPROTECT(1);`	405	`UNPROTECT(1);`
369	`PROTECT(OS->R_fcall = lang2(fn, R_NilValue));`	406	`PROTECT(OS->R_fcall = lang2(fn, R_NilValue));`
370	`PROTECT(par = coerceVector(par, REALSXP));`	407	`PROTECT(par = coerceVector(par, REALSXP));`
371	`if (!isNull(gr)) {`	408	`if (!isNull(gr)) {`
372	`if (!isFunction(gr)) error("gr is not a function");`	409	`if (!isFunction(gr)) error("gr is not a function");`
373	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`	410	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`
374	`} else {`	411	`} else {`
375	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`	412	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`
376	`}`	413	`}`
377	`ndeps = getListElement(options, "ndeps");`	414	`ndeps = getListElement(options, "ndeps");`
378	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`	415	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`
379	`OS->ndeps = vect(npar);`	416	`OS->ndeps = vect(npar);`
380	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`	417	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`
381	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`	418	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`
382	`UNPROTECT(1);`	419	`UNPROTECT(1);`
383	`PROTECT(ans = allocMatrix(REALSXP, npar, npar));`	420	`PROTECT(ans = allocMatrix(REALSXP, npar, npar));`
384	`dpar = vect(npar);`	421	`dpar = vect(npar);`
385	`for (i = 0; i < npar; i++)`	422	`for (i = 0; i < npar; i++)`
386	`dpar[i] = REAL(par)[i] / (OS->parscale[i]);`	423	`dpar[i] = REAL(par)[i] / (OS->parscale[i]);`
387	`df1 = vect(npar);`	424	`df1 = vect(npar);`
388	`df2 = vect(npar);`	425	`df2 = vect(npar);`
389	`for (i = 0; i < npar; i++) {`	426	`for (i = 0; i < npar; i++) {`
390	`eps = OS->ndeps[i]/(OS->parscale[i]);`	427	`eps = OS->ndeps[i]/(OS->parscale[i]);`
391	`dpar[i] = dpar[i] + eps;`	428	`dpar[i] = dpar[i] + eps;`
392	`fmingr(npar, dpar, df1, (void *)OS);`	429	`fmingr(npar, dpar, df1, (void *)OS);`
393	`dpar[i] = dpar[i] - 2 * eps;`	430	`dpar[i] = dpar[i] - 2 * eps;`
394	`fmingr(npar, dpar, df2, (void *)OS);`	431	`fmingr(npar, dpar, df2, (void *)OS);`
395	`for (j = 0; j < npar; j++)`	432	`for (j = 0; j < npar; j++)`
396	`REAL(ans)[i * npar + j] = (OS->fnscale) * (df1[j] - df2[j])/`	433	`REAL(ans)[i * npar + j] = (OS->fnscale) * (df1[j] - df2[j])/`
397	`(2 * eps * (OS->parscale[i]) * (OS->parscale[j]));`	434	`(2 * eps * (OS->parscale[i]) * (OS->parscale[j]));`
398	`dpar[i] = dpar[i] + eps;`	435	`dpar[i] = dpar[i] + eps;`
399	`}`	436	`}`
400	`vmaxset(vmax);`	437	`vmaxset(vmax);`
401	`UNPROTECT(4);`	438	`UNPROTECT(4);`
402	`return ans;`	439	`return ans;`
403	`}`	440	`}`
404		441
405		442
406	`static double ** matrix(int nrh, int nch)`	443	`static double ** matrix(int nrh, int nch)`
407	`{`	444	`{`
408	`int i;`	445	`int i;`
409	`double **m;`	446	`double **m;`
410		447
411	`m = (double *) R_alloc((nrh + 1), sizeof(double ));`	448	`m = (double *) R_alloc((nrh + 1), sizeof(double ));`
412	`for (i = 0; i <= nrh; i++)`	449	`for (i = 0; i <= nrh; i++)`
413	`m[i] = (double*) R_alloc((nch + 1), sizeof(double));`	450	`m[i] = (double*) R_alloc((nch + 1), sizeof(double));`
414	`return m;`	451	`return m;`
415	`}`	452	`}`
416		453
417	`static double ** Lmatrix(int n)`	454	`static double ** Lmatrix(int n)`
418	`{`	455	`{`
419	`int i;`	456	`int i;`
420	`double **m;`	457	`double **m;`
421		458
422	`m = (double *) R_alloc(n, sizeof(double ));`	459	`m = (double *) R_alloc(n, sizeof(double ));`
423	`for (i = 0; i < n; i++)`	460	`for (i = 0; i < n; i++)`
424	`m[i] = (double *) R_alloc((i + 1), sizeof(double));`	461	`m[i] = (double *) R_alloc((i + 1), sizeof(double));`
425	`return m;`	462	`return m;`
426	`}`	463	`}`
427		464
428		465
429		466
430	`#define stepredn 0.2`	467	`#define stepredn 0.2`
431	`#define acctol 0.0001`	468	`#define acctol 0.0001`
432	`#define reltest 10.0`	469	`#define reltest 10.0`
433		470
434		471
435	`/* BFGS variable-metric method, based on Pascal code`	472	`/* BFGS variable-metric method, based on Pascal code`
436	in J.C. Nash, `Compact Numerical Methods for Computers', 2nd edition,	473	in J.C. Nash, `Compact Numerical Methods for Computers', 2nd edition,
437	`converted by p2c then re-crafted by B.D. Ripley */`	474	`converted by p2c then re-crafted by B.D. Ripley */`
438		475
439	`void`	476	`void`
440	`vmmin(int n0, double b, double Fmin, optimfn fminfn, optimgr fmingr,`	477	`vmmin(int n0, double b, double Fmin, optimfn fminfn, optimgr fmingr,`
441	`int maxit, int trace, int *mask,`	478	`int maxit, int trace, int *mask,`
442	`double abstol, double reltol, int nREPORT, void *ex,`	479	`double abstol, double reltol, int nREPORT, void *ex,`
443	`int fncount, int grcount, int *fail)`	480	`int fncount, int grcount, int *fail)`
444	`{`	481	`{`
445	`Rboolean accpoint, enough;`	482	`Rboolean accpoint, enough;`
446	`double g, t, X, c, **B;`	483	`double g, t, X, c, **B;`
447	`int count, funcount, gradcount;`	484	`int count, funcount, gradcount;`
448	`double f, gradproj;`	485	`double f, gradproj;`
449	`int i, j, ilast, iter = 0;`	486	`int i, j, ilast, iter = 0;`
450	`double s, steplength;`	487	`double s, steplength;`
451	`double D1, D2;`	488	`double D1, D2;`
452	`int n, *l;`	489	`int n, *l;`
453		490
454	`if (maxit <= 0) {`	491	`if (maxit <= 0) {`
455	`*fail = 0;`	492	`*fail = 0;`
456	`*Fmin = fminfn(n0, b, ex);`	493	`*Fmin = fminfn(n0, b, ex);`
457	`fncount = grcount = 0;`	494	`fncount = grcount = 0;`
458	`return;`	495	`return;`
459	`}`	496	`}`
460		497
461	`if (nREPORT <= 0)`	498	`if (nREPORT <= 0)`
462	`error("REPORT must be > 0 (method = \"BFGS\")");`	499	`error("REPORT must be > 0 (method = \"BFGS\")");`
463	`l = (int *) R_alloc(n0, sizeof(int));`	500	`l = (int *) R_alloc(n0, sizeof(int));`

The R Project SVN R

(root)/trunk/src/main/optim.c @ 28362 – Rev 19233 → 21689