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

Rev	Author	Line No.	Line
7345	ripley	1	`/*`
		2	`* R : A Computer Language for Statistical Data Analysis`
18016	ripley	3	`* Copyright (C) 1999-2002 the R Development Core Team`
7345	ripley	4	`*`
		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`
		7	`* the Free Software Foundation; either version 2 of the License, or`
		8	`* (at your option) any later version.`
		9	`*`
		10	`* This program is distributed in the hope that it will be useful,`
		11	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`
		12	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`
		13	`* GNU General Public License for more details.`
		14	`*`
		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`
		17	`* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA`
		18	`*/`
		19
12778	pd	20	`#include <Defn.h>`
		21	`#include <R_ext/Random.h> /* for the random number generation in`
11497	hornik	22	`samin() */`
12778	pd	23	`#include <R_ext/Applic.h> /* setulb() */`
7345	ripley	24
		25	`static SEXP getListElement(SEXP list, char *str)`
		26	`{`
		27	`SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);`
		28	`int i;`
8186	ripley	29
7345	ripley	30	`for (i = 0; i < length(list); i++)`
10172	luke	31	`if (strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {`
		32	`elmt = VECTOR_ELT(list, i);`
7345	ripley	33	`break;`
		34	`}`
		35	`return elmt;`
		36	`}`
		37
		38	`static double * vect(int n)`
		39	`{`
		40	`return (double *)R_alloc(n, sizeof(double));`
		41	`}`
		42
8186	ripley	43	`typedef struct opt_struct`
7345	ripley	44	`{`
		45	`SEXP R_fcall; /* function */`
		46	`SEXP R_gcall; /* gradient */`
		47	`SEXP R_env; /* where to evaluate the calls */`
		48	`double* ndeps; /* tolerances for numerical derivatives */`
		49	`double fnscale; /* scaling for objective */`
		50	`double* parscale;/* scaling for parameters */`
7974	ripley	51	`int usebounds;`
		52	`double* lower, *upper;`
7345	ripley	53	`} opt_struct, *OptStruct;`
		54
		55
7464	ripley	56
18016	ripley	57	`static double fminfn(int n, double p, void ex)`
7345	ripley	58	`{`
		59	`SEXP s, x;`
		60	`int i;`
		61	`double val;`
18016	ripley	62	`OptStruct OS = (OptStruct) ex;`
10853	ripley	63	`PROTECT_INDEX ipx;`
8186	ripley	64
7345	ripley	65	`PROTECT(x = allocVector(REALSXP, n));`
8316	ripley	66	`for (i = 0; i < n; i++) {`
		67	`if (!R_FINITE(p[i])) error("non-finite value supplied by optim");`
		68	`REAL(x)[i] = p[i] * (OS->parscale[i]);`
		69	`}`
10172	luke	70	`SETCADR(OS->R_fcall, x);`
10853	ripley	71	`PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);`
		72	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
7345	ripley	73	`val = REAL(s)[0]/(OS->fnscale);`
		74	`UNPROTECT(2);`
		75	`return val;`
		76	`}`
		77
18016	ripley	78	`static void fmingr(int n, double p, double df, void *ex)`
7345	ripley	79	`{`
		80	`SEXP s, x;`
		81	`int i;`
7974	ripley	82	`double val1, val2, eps, epsused, tmp;`
18016	ripley	83	`OptStruct OS = (OptStruct) ex;`
10853	ripley	84	`PROTECT_INDEX ipx;`
7345	ripley	85
7510	ripley	86	`if (!isNull(OS->R_gcall)) { /* analytical derivatives */`
7345	ripley	87	`PROTECT(x = allocVector(REALSXP, n));`
8316	ripley	88	`for (i = 0; i < n; i++) {`
10196	maechler	89	`if (!R_FINITE(p[i])) error("non-finite value supplied by optim");`
7345	ripley	90	`REAL(x)[i] = p[i] * (OS->parscale[i]);`
8316	ripley	91	`}`
10172	luke	92	`SETCADR(OS->R_gcall, x);`
10853	ripley	93	`PROTECT_WITH_INDEX(s = eval(OS->R_gcall, OS->R_env), &ipx);`
		94	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
19233	maechler	95	`if(LENGTH(s) != n)`
		96	`error("gradient in optim evaluated to length %d not %d",`
11347	ripley	97	`LENGTH(s), n);`
8186	ripley	98	`for (i = 0; i < n; i++)`
7345	ripley	99	`df[i] = REAL(s)[i] * (OS->parscale[i])/(OS->fnscale);`
		100	`UNPROTECT(2);`
		101	`} else { /* numerical derivatives */`
		102	`PROTECT(x = allocVector(REALSXP, n));`
		103	`for (i = 0; i < n; i++) REAL(x)[i] = p[i] * (OS->parscale[i]);`
10172	luke	104	`SETCADR(OS->R_fcall, x);`
7974	ripley	105	`if(OS->usebounds == 0) {`
		106	`for (i = 0; i < n; i++) {`
		107	`eps = OS->ndeps[i];`
		108	`REAL(x)[i] = (p[i] + eps) * (OS->parscale[i]);`
10172	luke	109	`SETCADR(OS->R_fcall, x);`
10853	ripley	110	`PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);`
		111	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
7974	ripley	112	`val1 = REAL(s)[0]/(OS->fnscale);`
		113	`REAL(x)[i] = (p[i] - eps) * (OS->parscale[i]);`
10172	luke	114	`SETCADR(OS->R_fcall, x);`
10853	ripley	115	`REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);`
		116	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
7974	ripley	117	`val2 = REAL(s)[0]/(OS->fnscale);`
		118	`df[i] = (val1 - val2)/(2 * eps);`
10196	maechler	119	`#define DO_df_x \`
		120	`if(!R_FINITE(df[i])) \`
		121	`error("non-finite finite-difference value [%d]", i);\`
		122	`REAL(x)[i] = p[i] * (OS->parscale[i])`
		123
		124	`DO_df_x;`
10853	ripley	125	`UNPROTECT(1);`
7974	ripley	126	`}`
10196	maechler	127	`} else { /* usebounds */`
7974	ripley	128	`for (i = 0; i < n; i++) {`
		129	`epsused = eps = OS->ndeps[i];`
10853	ripley	130	`tmp = p[i] + eps;`
7974	ripley	131	`if (tmp > OS->upper[i]) {`
		132	`tmp = OS->upper[i];`
		133	`epsused = tmp - p[i] ;`
		134	`}`
		135	`REAL(x)[i] = tmp * (OS->parscale[i]);`
10172	luke	136	`SETCADR(OS->R_fcall, x);`
10853	ripley	137	`PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);`
		138	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
7974	ripley	139	`val1 = REAL(s)[0]/(OS->fnscale);`
		140	`tmp = p[i] - eps;`
		141	`if (tmp < OS->lower[i]) {`
		142	`tmp = OS->lower[i];`
		143	`eps = p[i] - tmp;`
		144	`}`
		145	`REAL(x)[i] = tmp * (OS->parscale[i]);`
10172	luke	146	`SETCADR(OS->R_fcall, x);`
10853	ripley	147	`REPROTECT(s = eval(OS->R_fcall, OS->R_env), ipx);`
		148	`REPROTECT(s = coerceVector(s, REALSXP), ipx);`
7974	ripley	149	`val2 = REAL(s)[0]/(OS->fnscale);`
		150	`df[i] = (val1 - val2)/(epsused + eps);`
10196	maechler	151
		152	`DO_df_x;`
10853	ripley	153	`UNPROTECT(1);`
7974	ripley	154	`}`
7345	ripley	155	`}`
		156	`UNPROTECT(1); /* x */`
		157	`}`
		158	`}`
		159
21689	ripley	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
7345	ripley	190	`/* par fn gr method options */`
		191	`SEXP do_optim(SEXP call, SEXP op, SEXP args, SEXP rho)`
		192	`{`
7464	ripley	193	`SEXP par, fn, gr, method, options, tmp, slower, supper;`
7345	ripley	194	`SEXP res, value, counts, conv;`
8186	ripley	195	`int i, npar=0, *mask, trace, maxit, fncount, grcount, nREPORT, tmax;`
7345	ripley	196	`int ifail = 0;`
8186	ripley	197	`double dpar, opar, val, abstol, reltol, temp;`
7345	ripley	198	`char *tn;`
		199	`OptStruct OS;`
7975	ripley	200	`char *vmax;`
7345	ripley	201
		202	`checkArity(op, args);`
7975	ripley	203	`vmax = vmaxget();`
7345	ripley	204	`OS = (OptStruct) R_alloc(1, sizeof(opt_struct));`
7974	ripley	205	`OS->usebounds = 0;`
7345	ripley	206	`OS->R_env = rho;`
		207	`par = CAR(args);`
		208	`args = CDR(args); fn = CAR(args);`
7510	ripley	209	`if (!isFunction(fn)) errorcall(call, "fn is not a function");`
7345	ripley	210	`args = CDR(args); gr = CAR(args);`
		211	`args = CDR(args); method = CAR(args);`
8186	ripley	212	`if (!isString(method)\|\| LENGTH(method) != 1)`
7345	ripley	213	`errorcall(call, "invalid method argument");`
10172	luke	214	`tn = CHAR(STRING_ELT(method, 0));`
7345	ripley	215	`args = CDR(args); options = CAR(args);`
		216	`PROTECT(OS->R_fcall = lang2(fn, R_NilValue));`
7974	ripley	217	`PROTECT(par = coerceVector(duplicate(par), REALSXP));`
7345	ripley	218	`npar = LENGTH(par);`
		219	`dpar = vect(npar);`
		220	`opar = vect(npar);`
		221	`trace = asInteger(getListElement(options, "trace"));`
		222	`OS->fnscale = asReal(getListElement(options, "fnscale"));`
		223	`tmp = getListElement(options, "parscale");`
8186	ripley	224	`if (LENGTH(tmp) != npar)`
7345	ripley	225	`errorcall(call, "parscale is of the wrong length");`
		226	`PROTECT(tmp = coerceVector(tmp, REALSXP));`
		227	`OS->parscale = vect(npar);`
		228	`for (i = 0; i < npar; i++) OS->parscale[i] = REAL(tmp)[i];`
		229	`UNPROTECT(1);`
8186	ripley	230	`for (i = 0; i < npar; i++)`
7345	ripley	231	`dpar[i] = REAL(par)[i] / (OS->parscale[i]);`
7464	ripley	232	`PROTECT(res = allocVector(VECSXP, 5));`
7345	ripley	233	`PROTECT(value = allocVector(REALSXP, 1));`
		234	`PROTECT(counts = allocVector(INTSXP, 2));`
		235	`PROTECT(conv = allocVector(INTSXP, 1));`
		236	`abstol = asReal(getListElement(options, "abstol"));`
8186	ripley	237	`reltol = asReal(getListElement(options, "reltol"));`
		238	`maxit = asInteger(getListElement(options, "maxit"));`
7510	ripley	239	`if (maxit == NA_INTEGER) error("maxit is not an integer");`
7345	ripley	240
		241	`if (strcmp(tn, "Nelder-Mead") == 0) {`
10887	maechler	242	`double alpha, beta, gamm;`
7345	ripley	243
		244	`alpha = asReal(getListElement(options, "alpha"));`
		245	`beta = asReal(getListElement(options, "beta"));`
10887	maechler	246	`gamm = asReal(getListElement(options, "gamma"));`
19233	maechler	247	`nmmin(npar, dpar, opar, &val, fminfn, &ifail, abstol, reltol,`
18016	ripley	248	`(void *)OS, alpha, beta, gamm, trace, &fncount, maxit);`
8186	ripley	249	`for (i = 0; i < npar; i++)`
7345	ripley	250	`REAL(par)[i] = opar[i] * (OS->parscale[i]);`
		251	`grcount = NA_INTEGER;`
8186	ripley	252
		253	`}`
21689	ripley	254	`else if (strcmp(tn, "SANN") == 0) {`
		255	`tmax = asInteger(getListElement(options, "tmax"));`
		256	`temp = asReal(getListElement(options, "temp"));`
		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	`}`
		264	`samin (npar, dpar, &val, fminfn, maxit, tmax, temp, trace, (void *)OS);`
		265	`for (i = 0; i < npar; i++)`
		266	`REAL(par)[i] = dpar[i] * (OS->parscale[i]);`
		267	`fncount = maxit;`
		268	`grcount = NA_INTEGER;`
		269	`UNPROTECT(1); /* OS->R_gcall */`
8186	ripley	270
7345	ripley	271	`} else if (strcmp(tn, "BFGS") == 0) {`
		272	`SEXP ndeps;`
8186	ripley	273
7345	ripley	274	`nREPORT = asInteger(getListElement(options, "REPORT"));`
7510	ripley	275	`if (!isNull(gr)) {`
		276	`if (!isFunction(gr)) error("gr is not a function");`
7345	ripley	277	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`
		278	`} else {`
		279	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`
8186	ripley	280	`ndeps = getListElement(options, "ndeps");`
7510	ripley	281	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`
7345	ripley	282	`OS->ndeps = vect(npar);`
		283	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`
		284	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`
		285	`UNPROTECT(1);`
		286	`}`
		287	`mask = (int *) R_alloc(npar, sizeof(int));`
		288	`for (i = 0; i < npar; i++) mask[i] = 1;`
19233	maechler	289	`vmmin(npar, dpar, &val, fminfn, fmingr, maxit, trace, mask, abstol,`
18016	ripley	290	`reltol, nREPORT, (void *)OS, &fncount, &grcount, &ifail);`
8186	ripley	291	`for (i = 0; i < npar; i++)`
7345	ripley	292	`REAL(par)[i] = dpar[i] * (OS->parscale[i]);`
		293	`UNPROTECT(1); /* OS->R_gcall */`
		294	`} else if (strcmp(tn, "CG") == 0) {`
		295	`int type;`
		296	`SEXP ndeps;`
		297
		298	`type = asInteger(getListElement(options, "type"));`
7510	ripley	299	`if (!isNull(gr)) {`
		300	`if (!isFunction(gr)) error("gr is not a function");`
7345	ripley	301	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`
		302	`} else {`
		303	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`
8186	ripley	304	`ndeps = getListElement(options, "ndeps");`
7510	ripley	305	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`
7345	ripley	306	`OS->ndeps = vect(npar);`
		307	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`
		308	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`
		309	`UNPROTECT(1);`
		310	`}`
19233	maechler	311	`cgmin(npar, dpar, opar, &val, fminfn, fmingr, &ifail, abstol,`
18016	ripley	312	`reltol, (void *)OS, type, trace, &fncount, &grcount, maxit);`
8186	ripley	313	`for (i = 0; i < npar; i++)`
7345	ripley	314	`REAL(par)[i] = opar[i] * (OS->parscale[i]);`
		315	`UNPROTECT(1); /* OS->R_gcall */`
10853	ripley	316
7464	ripley	317	`} else if (strcmp(tn, "L-BFGS-B") == 0) {`
		318	`SEXP ndeps, smsg;`
		319	`double lower = vect(npar), upper = vect(npar);`
		320	`int lmm, nbd = (int ) R_alloc(npar, sizeof(int));`
		321	`double factr, pgtol;`
		322	`char msg[60];`
		323
10853	ripley	324	`nREPORT = asInteger(getListElement(options, "REPORT"));`
7464	ripley	325	`factr = asReal(getListElement(options, "factr"));`
		326	`pgtol = asReal(getListElement(options, "pgtol"));`
		327	`lmm = asInteger(getListElement(options, "lmm"));`
7510	ripley	328	`if (!isNull(gr)) {`
		329	`if (!isFunction(gr)) error("gr is not a function");`
7464	ripley	330	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`
		331	`} else {`
		332	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`
8186	ripley	333	`ndeps = getListElement(options, "ndeps");`
7510	ripley	334	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`
7464	ripley	335	`OS->ndeps = vect(npar);`
		336	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`
		337	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`
		338	`UNPROTECT(1);`
		339	`}`
		340	`args = CDR(args); slower = CAR(args); /* coerce in calling code */`
		341	`args = CDR(args); supper = CAR(args);`
		342	`for (i = 0; i < npar; i++) {`
7972	ripley	343	`lower[i] = REAL(slower)[i] / (OS->parscale[i]);`
		344	`upper[i] = REAL(supper)[i] / (OS->parscale[i]);`
7464	ripley	345	`if (!R_FINITE(lower[i])) {`
		346	`if (!R_FINITE(upper[i])) nbd[i] = 0; else nbd[i] = 3;`
		347	`} else {`
		348	`if (!R_FINITE(upper[i])) nbd[i] = 1; else nbd[i] = 2;`
		349	`}`
		350	`}`
7974	ripley	351	`OS->usebounds = 1;`
		352	`OS->lower = lower;`
		353	`OS->upper = upper;`
19233	maechler	354	`lbfgsb(npar, lmm, dpar, lower, upper, nbd, &val, fminfn, fmingr,`
		355	`&ifail, (void *)OS, factr, pgtol, &fncount, &grcount,`
18016	ripley	356	`maxit, msg, trace, nREPORT);`
8186	ripley	357	`for (i = 0; i < npar; i++)`
7464	ripley	358	`REAL(par)[i] = dpar[i] * (OS->parscale[i]);`
		359	`UNPROTECT(1); /* OS->R_gcall */`
		360	`PROTECT(smsg = allocVector(STRSXP, 1));`
19105	ripley	361	`SET_STRING_ELT(smsg, 0, mkChar(msg));`
10172	luke	362	`SET_VECTOR_ELT(res, 4, smsg);`
7464	ripley	363	`UNPROTECT(1);`
8186	ripley	364	`} else`
7464	ripley	365	`errorcall(call, "unknown method");`
		366
7345	ripley	367	`REAL(value)[0] = val * (OS->fnscale);`
10172	luke	368	`SET_VECTOR_ELT(res, 0, par); SET_VECTOR_ELT(res, 1, value);`
7345	ripley	369	`INTEGER(counts)[0] = fncount; INTEGER(counts)[1] = grcount;`
10172	luke	370	`SET_VECTOR_ELT(res, 2, counts);`
7345	ripley	371	`INTEGER(conv)[0] = ifail;`
10172	luke	372	`SET_VECTOR_ELT(res, 3, conv);`
7975	ripley	373	`vmaxset(vmax);`
7345	ripley	374	`UNPROTECT(6);`
		375	`return res;`
		376	`}`
		377
		378	`/* par fn gr options */`
		379	`SEXP do_optimhess(SEXP call, SEXP op, SEXP args, SEXP rho)`
		380	`{`
		381	`SEXP par, fn, gr, options, tmp, ndeps, ans;`
		382	`OptStruct OS;`
		383	`int npar, i , j;`
		384	`double dpar, df1, *df2, eps;`
7975	ripley	385	`char *vmax;`
7345	ripley	386
		387	`checkArity(op, args);`
7975	ripley	388	`vmax = vmaxget();`
7345	ripley	389	`OS = (OptStruct) R_alloc(1, sizeof(opt_struct));`
8988	luke	390	`OS->usebounds = 0;`
7345	ripley	391	`OS->R_env = rho;`
		392	`par = CAR(args);`
		393	`npar = LENGTH(par);`
		394	`args = CDR(args); fn = CAR(args);`
7510	ripley	395	`if (!isFunction(fn)) errorcall(call, "fn is not a function");`
7345	ripley	396	`args = CDR(args); gr = CAR(args);`
		397	`args = CDR(args); options = CAR(args);`
		398	`OS->fnscale = asReal(getListElement(options, "fnscale"));`
		399	`tmp = getListElement(options, "parscale");`
8186	ripley	400	`if (LENGTH(tmp) != npar)`
7345	ripley	401	`errorcall(call, "parscale is of the wrong length");`
		402	`PROTECT(tmp = coerceVector(tmp, REALSXP));`
		403	`OS->parscale = vect(npar);`
		404	`for (i = 0; i < npar; i++) OS->parscale[i] = REAL(tmp)[i];`
		405	`UNPROTECT(1);`
		406	`PROTECT(OS->R_fcall = lang2(fn, R_NilValue));`
		407	`PROTECT(par = coerceVector(par, REALSXP));`
7510	ripley	408	`if (!isNull(gr)) {`
		409	`if (!isFunction(gr)) error("gr is not a function");`
7345	ripley	410	`PROTECT(OS->R_gcall = lang2(gr, R_NilValue));`
		411	`} else {`
		412	`PROTECT(OS->R_gcall = R_NilValue); /* for balance */`
		413	`}`
8186	ripley	414	`ndeps = getListElement(options, "ndeps");`
7510	ripley	415	`if (LENGTH(ndeps) != npar) error("ndeps is of the wrong length");`
7345	ripley	416	`OS->ndeps = vect(npar);`
		417	`PROTECT(ndeps = coerceVector(ndeps, REALSXP));`
		418	`for (i = 0; i < npar; i++) OS->ndeps[i] = REAL(ndeps)[i];`
		419	`UNPROTECT(1);`
		420	`PROTECT(ans = allocMatrix(REALSXP, npar, npar));`
		421	`dpar = vect(npar);`
8186	ripley	422	`for (i = 0; i < npar; i++)`
7345	ripley	423	`dpar[i] = REAL(par)[i] / (OS->parscale[i]);`
		424	`df1 = vect(npar);`
		425	`df2 = vect(npar);`
		426	`for (i = 0; i < npar; i++) {`
		427	`eps = OS->ndeps[i]/(OS->parscale[i]);`
		428	`dpar[i] = dpar[i] + eps;`
18016	ripley	429	`fmingr(npar, dpar, df1, (void *)OS);`
7345	ripley	430	`dpar[i] = dpar[i] - 2 * eps;`
18016	ripley	431	`fmingr(npar, dpar, df2, (void *)OS);`
7345	ripley	432	`for (j = 0; j < npar; j++)`
		433	`REAL(ans)[i * npar + j] = (OS->fnscale) * (df1[j] - df2[j])/`
		434	`(2 * eps * (OS->parscale[i]) * (OS->parscale[j]));`
		435	`dpar[i] = dpar[i] + eps;`
		436	`}`
7975	ripley	437	`vmaxset(vmax);`
7345	ripley	438	`UNPROTECT(4);`
		439	`return ans;`
		440	`}`
		441
		442
		443	`static double ** matrix(int nrh, int nch)`
		444	`{`
		445	`int i;`
		446	`double **m;`
		447
		448	`m = (double *) R_alloc((nrh + 1), sizeof(double ));`
8186	ripley	449	`for (i = 0; i <= nrh; i++)`
7345	ripley	450	`m[i] = (double*) R_alloc((nch + 1), sizeof(double));`
		451	`return m;`
		452	`}`
		453
		454	`static double ** Lmatrix(int n)`
		455	`{`
		456	`int i;`
		457	`double **m;`
		458
		459	`m = (double *) R_alloc(n, sizeof(double ));`
8186	ripley	460	`for (i = 0; i < n; i++)`
7345	ripley	461	`m[i] = (double *) R_alloc((i + 1), sizeof(double));`
		462	`return m;`
		463	`}`
		464
		465
		466
		467	`#define stepredn 0.2`
		468	`#define acctol 0.0001`
		469	`#define reltest 10.0`
		470
		471
		472	`/* BFGS variable-metric method, based on Pascal code`
		473	in J.C. Nash, `Compact Numerical Methods for Computers', 2nd edition,
		474	`converted by p2c then re-crafted by B.D. Ripley */`
		475
18031	ripley	476	`void`
19233	maechler	477	`vmmin(int n0, double b, double Fmin, optimfn fminfn, optimgr fmingr,`
18016	ripley	478	`int maxit, int trace, int *mask,`
		479	`double abstol, double reltol, int nREPORT, void *ex,`
7345	ripley	480	`int fncount, int grcount, int *fail)`
		481	`{`
10470	maechler	482	`Rboolean accpoint, enough;`
7345	ripley	483	`double g, t, X, c, **B;`
		484	`int count, funcount, gradcount;`
		485	`double f, gradproj;`
		486	`int i, j, ilast, iter = 0;`
		487	`double s, steplength;`
		488	`double D1, D2;`
		489	`int n, *l;`
		490
18222	ripley	491	`if (maxit <= 0) {`
18295	ripley	492	`*fail = 0;`
18222	ripley	493	`*Fmin = fminfn(n0, b, ex);`
		494	`fncount = grcount = 0;`
		495	`return;`
		496	`}`
		497
10196	maechler	498	`if (nREPORT <= 0)`
		499	`error("REPORT must be > 0 (method = \"BFGS\")");`
7345	ripley	500	`l = (int *) R_alloc(n0, sizeof(int));`
		501	`n = 0;`
		502	`for (i = 0; i < n0; i++) if (mask[i]) l[n++] = i;`
		503	`g = vect(n0);`
		504	`t = vect(n);`
		505	`X = vect(n);`
		506	`c = vect(n);`
		507	`B = Lmatrix(n);`
18222	ripley	508	`f = fminfn(n0, b, ex);`
10196	maechler	509	`if (!R_FINITE(f))`
		510	`error("initial value in vmmin is not finite");`
7345	ripley	511	`if (trace) Rprintf("initial value %f \n", f);`
		512	`*Fmin = f;`
		513	`funcount = gradcount = 1;`
18222	ripley	514	`fmingr(n0, b, g, ex);`
7345	ripley	515	`iter++;`
		516	`ilast = gradcount;`
8186	ripley	517
7345	ripley	518	`do {`
		519	`if (ilast == gradcount) {`
		520	`for (i = 0; i < n; i++) {`
		521	`for (j = 0; j < i; j++) B[i][j] = 0.0;`
		522	`B[i][i] = 1.0;`
		523	`}`
		524	`}`
		525	`for (i = 0; i < n; i++) {`
		526	`X[i] = b[l[i]];`
		527	`c[i] = g[l[i]];`
		528	`}`
		529	`gradproj = 0.0;`
		530	`for (i = 0; i < n; i++) {`
		531	`s = 0.0;`
		532	`for (j = 0; j <= i; j++) s -= B[i][j] * g[l[j]];`
		533	`for (j = i + 1; j < n; j++) s -= B[j][i] * g[l[j]];`
		534	`t[i] = s;`
		535	`gradproj += s * g[l[i]];`
		536	`}`
		537
		538	`if (gradproj < 0.0) { /* search direction is downhill */`
		539	`steplength = 1.0;`
10470	maechler	540	`accpoint = FALSE;`
7345	ripley	541	`do {`
		542	`count = 0;`
		543	`for (i = 0; i < n; i++) {`
		544	`b[l[i]] = X[i] + steplength * t[i];`
		545	`if (reltest + X[i] == reltest + b[l[i]]) /* no change */`
		546	`count++;`
		547	`}`
		548	`if (count < n) {`
18222	ripley	549	`f = fminfn(n0, b, ex);`
7345	ripley	550	`funcount++;`
8186	ripley	551	`accpoint = R_FINITE(f) &&`
7345	ripley	552	`(f <= Fmin + gradproj steplength * acctol);`
		553	`if (!accpoint) {`
		554	`steplength *= stepredn;`
		555	`}`
		556	`}`
		557	`} while (!(count == n \|\| accpoint));`
8186	ripley	558	`enough = (f > abstol) &&`
7345	ripley	559	`fabs(f - Fmin) > reltol (fabs(*Fmin) + reltol);`
		560	`/* stop if value if small or if relative change is low */`
		561	`if (!enough) {`
		562	`count = n;`
		563	`*Fmin = f;`
		564	`}`
		565	`if (count < n) {/* making progress */`
		566	`*Fmin = f;`
18222	ripley	567	`fmingr(n0, b, g, ex);`
7345	ripley	568	`gradcount++;`
		569	`iter++;`
		570	`D1 = 0.0;`
		571	`for (i = 0; i < n; i++) {`
		572	`t[i] = steplength * t[i];`
		573	`c[i] = g[l[i]] - c[i];`
		574	`D1 += t[i] * c[i];`
		575	`}`
		576	`if (D1 > 0) {`
		577	`D2 = 0.0;`
		578	`for (i = 0; i < n; i++) {`
		579	`s = 0.0;`
		580	`for (j = 0; j <= i; j++)`
		581	`s += B[i][j] * c[j];`
		582	`for (j = i + 1; j < n; j++)`
		583	`s += B[j][i] * c[j];`
		584	`X[i] = s;`
		585	`D2 += s * c[i];`
		586	`}`
		587	`D2 = 1.0 + D2 / D1;`
		588	`for (i = 0; i < n; i++) {`
		589	`for (j = 0; j <= i; j++)`
10196	maechler	590	`B[i][j] += (D2 * t[i] * t[j]`
		591	`- X[i] * t[j] - t[i] * X[j]) / D1;`
7345	ripley	592	`}`
		593	`} else { /* D1 < 0 */`
		594	`ilast = gradcount;`
		595	`}`
		596	`} else { /* no progress */`
		597	`if (ilast < gradcount) {`
		598	`count = 0;`
		599	`ilast = gradcount;`
		600	`}`
		601	`}`
		602	`} else { /* uphill search */`
		603	`count = 0;`
		604	`if (ilast == gradcount) count = n;`
		605	`else ilast = gradcount;`
		606	`/* Resets unless has just been reset */`
		607	`}`
10196	maechler	608	`if (trace && (iter % nREPORT == 0))`
7345	ripley	609	`Rprintf("iter%4d value %f\n", iter, f);`
		610	`if (iter >= maxit) break;`
		611	`if (gradcount - ilast > 2 * n)`
		612	`ilast = gradcount; /* periodic restart */`
		613	`} while (count != n \|\| ilast != gradcount);`
		614	`if (trace) {`
		615	`Rprintf("final value %f \n", *Fmin);`
		616	`if (iter < maxit) Rprintf("converged\n");`
		617	`else Rprintf("stopped after %i iterations\n", iter);`
		618	`}`
18295	ripley	619	`*fail = (iter < maxit) ? 0 : 1;`
7345	ripley	620	`*fncount = funcount;`
		621	`*grcount = gradcount;`
		622	`}`
		623
		624
		625	`#define big 1.0e+35 /a very large number/`
		626
		627
10196	maechler	628	`/* Nelder-Mead */`
18016	ripley	629	`void nmmin(int n, double Bvec, double X, double *Fmin, optimfn fminfn,`
		630	`int fail, double abstol, double intol, void ex,`
19233	maechler	631	`double alpha, double bet, double gamm, int trace,`
7345	ripley	632	`int *fncount, int maxit)`
		633	`{`
		634	`char action[50];`
		635	`int C;`
10470	maechler	636	`Rboolean calcvert, shrinkfail = FALSE;`
7345	ripley	637	`double convtol, f;`
		638	`int funcount=0, H, i, j, L=0;`
		639	`int n1=0;`
		640	`double oldsize;`
		641	`double **P;`
		642	`double size, step, temp, trystep;`
		643	`char tstr[6];`
12976	pd	644	`double VH, VL, VR;`
7345	ripley	645
18222	ripley	646	`if (maxit <= 0) {`
		647	`*Fmin = fminfn(n, Bvec, ex);`
		648	`*fncount = 0;`
18295	ripley	649	`*fail = 0;`
18222	ripley	650	`return;`
		651	`}`
7510	ripley	652	`if (trace)`
7345	ripley	653	`Rprintf(" Nelder-Mead direct search function minimizer\n");`
		654	`P = matrix(n, n+1);`
10470	maechler	655	`*fail = FALSE;`
18016	ripley	656	`f = fminfn(n, Bvec, ex);`
7345	ripley	657	`if (!R_FINITE(f)) {`
		658	`error("Function cannot be evaluated at initial parameters");`
10470	maechler	659	`*fail = TRUE;`
7345	ripley	660	`} else {`
7510	ripley	661	`if (trace) Rprintf("Function value for initial parameters = %f\n", f);`
7345	ripley	662	`funcount = 1;`
		663	`convtol = intol * (fabs(f) + intol);`
7510	ripley	664	`if (trace) Rprintf(" Scaled convergence tolerance is %g\n", convtol);`
7345	ripley	665	`n1 = n + 1;`
		666	`C = n + 2;`
		667	`P[n1 - 1][0] = f;`
		668	`for (i = 0; i < n; i++)`
		669	`P[i][0] = Bvec[i];`
		670
		671	`L = 1;`
		672	`size = 0.0;`
		673
		674	`step = 0.0;`
		675	`for (i = 0; i < n; i++) {`
		676	`if (0.1 * fabs(Bvec[i]) > step)`
		677	`step = 0.1 * fabs(Bvec[i]);`
		678	`}`
7510	ripley	679	`if (step == 0.0) step = 0.1;`
		680	`if (trace) Rprintf("Stepsize computed as %f\n", step);`
7345	ripley	681	`for (j = 2; j <= n1; j++) {`
		682	`strcpy(action, "BUILD ");`
		683	`for (i = 0; i < n; i++)`
		684	`P[i][j - 1] = Bvec[i];`
		685
		686	`trystep = step;`
		687	`while (P[j - 2][j - 1] == Bvec[j - 2]) {`
		688	`P[j - 2][j - 1] = Bvec[j - 2] + trystep;`
		689	`trystep *= 10;`
		690	`}`
		691	`size += trystep;`
		692	`}`
		693	`oldsize = size;`
10470	maechler	694	`calcvert = TRUE;`
		695	`shrinkfail = FALSE;`
7345	ripley	696	`do {`
		697	`if (calcvert) {`
		698	`for (j = 0; j < n1; j++) {`
		699	`if (j + 1 != L) {`
		700	`for (i = 0; i < n; i++)`
		701	`Bvec[i] = P[i][j];`
18016	ripley	702	`f = fminfn(n, Bvec, ex);`
10217	pd	703	`if (!R_FINITE(f)) f = big;`
7345	ripley	704	`funcount++;`
		705	`P[n1 - 1][j] = f;`
		706	`}`
		707	`}`
10470	maechler	708	`calcvert = FALSE;`
7345	ripley	709	`}`
		710
		711	`VL = P[n1 - 1][L - 1];`
		712	`VH = VL;`
		713	`H = L;`
		714
		715	`for (j = 1; j <= n1; j++) {`
		716	`if (j != L) {`
		717	`f = P[n1 - 1][j - 1];`
		718	`if (f < VL) {`
		719	`L = j;`
		720	`VL = f;`
		721	`}`
		722	`if (f > VH) {`
		723	`H = j;`
		724	`VH = f;`
		725	`}`
		726	`}`
		727	`}`
		728
		729	`if (VH > VL + convtol && VL > abstol) {`
		730	`sprintf(tstr, "%5d", funcount);`
7510	ripley	731	`if (trace) Rprintf("%s%s %f %f\n", action, tstr, VH, VL);`
7345	ripley	732
		733	`for (i = 0; i < n; i++) {`
		734	`temp = -P[i][H - 1];`
		735	`for (j = 0; j < n1; j++)`
		736	`temp += P[i][j];`
		737	`P[i][C - 1] = temp / n;`
		738	`}`
		739	`for (i = 0; i < n; i++)`
		740	`Bvec[i] = (1.0 + alpha) * P[i][C - 1] - alpha * P[i][H - 1];`
18016	ripley	741	`f = fminfn(n, Bvec, ex);`
7345	ripley	742	`if (!R_FINITE(f)) f = big;`
		743	`funcount++;`
		744	`strcpy(action, "REFLECTION ");`
		745	`VR = f;`
		746	`if (VR < VL) {`
		747	`P[n1 - 1][C - 1] = f;`
		748	`for (i = 0; i < n; i++) {`
10887	maechler	749	`f = gamm * Bvec[i] + (1 - gamm) * P[i][C - 1];`
7345	ripley	750	`P[i][C - 1] = Bvec[i];`
		751	`Bvec[i] = f;`
		752	`}`
18016	ripley	753	`f = fminfn(n, Bvec, ex);`
7345	ripley	754	`if (!R_FINITE(f)) f = big;`
		755	`funcount++;`
		756	`if (f < VR) {`
		757	`for (i = 0; i < n; i++)`
		758	`P[i][H - 1] = Bvec[i];`
		759	`P[n1 - 1][H - 1] = f;`
		760	`strcpy(action, "EXTENSION ");`
		761	`} else {`
		762	`for (i = 0; i < n; i++)`
		763	`P[i][H - 1] = P[i][C - 1];`
		764	`P[n1 - 1][H - 1] = VR;`
		765	`}`
		766	`} else {`
10196	maechler	767	`strcpy(action, "HI-REDUCTION ");`
7345	ripley	768	`if (VR < VH) {`
		769	`for (i = 0; i < n; i++)`
		770	`P[i][H - 1] = Bvec[i];`
		771	`P[n1 - 1][H - 1] = VR;`
10196	maechler	772	`strcpy(action, "LO-REDUCTION ");`
7345	ripley	773	`}`
		774
		775	`for (i = 0; i < n; i++)`
19233	maechler	776	`Bvec[i] = (1 - bet) * P[i][H - 1] + bet * P[i][C - 1];`
18016	ripley	777	`f = fminfn(n, Bvec, ex);`
7345	ripley	778	`if (!R_FINITE(f)) f = big;`
		779	`funcount++;`
		780
		781	`if (f < P[n1 - 1][H - 1]) {`
		782	`for (i = 0; i < n; i++)`
		783	`P[i][H - 1] = Bvec[i];`
		784	`P[n1 - 1][H - 1] = f;`
		785	`} else {`
		786	`if (VR >= VH) {`
		787	`strcpy(action, "SHRINK ");`
10470	maechler	788	`calcvert = TRUE;`
7345	ripley	789	`size = 0.0;`
		790	`for (j = 0; j < n1; j++) {`
		791	`if (j + 1 != L) {`
		792	`for (i = 0; i < n; i++) {`
19233	maechler	793	`P[i][j] = bet * (P[i][j] - P[i][L - 1])`
		794	`+ P[i][L - 1];`
7345	ripley	795	`size += fabs(P[i][j] - P[i][L - 1]);`
		796	`}`
		797	`}`
		798	`}`
		799	`if (size < oldsize) {`
10470	maechler	800	`shrinkfail = FALSE;`
7345	ripley	801	`oldsize = size;`
		802	`} else {`
7510	ripley	803	`if (trace)`
7345	ripley	804	`Rprintf("Polytope size measure not decreased in shrink\n");`
10470	maechler	805	`shrinkfail = TRUE;`
7345	ripley	806	`}`
		807	`}`
		808	`}`
		809	`}`
		810	`}`
		811
8186	ripley	812	`} while (!(VH <= VL + convtol \|\| VL <= abstol \|\|`
7345	ripley	813	`shrinkfail \|\| funcount > maxit));`
		814
		815	`}`
		816
7510	ripley	817	`if (trace) {`
7345	ripley	818	`Rprintf("Exiting from Nelder Mead minimizer\n");`
		819	`Rprintf(" %d function evaluations used\n", funcount);`
		820	`}`
		821	`*Fmin = P[n1 - 1][L - 1];`
		822	`for (i = 0; i < n; i++) X[i] = P[i][L - 1];`
		823	`if (shrinkfail) *fail = 10;`
		824	`if (funcount > maxit) *fail = 1;`
		825	`*fncount = funcount;`
		826	`}`
		827
19233	maechler	828	`void cgmin(int n, double Bvec, double X, double *Fmin,`
18016	ripley	829	`optimfn fminfn, optimgr fmingr, int *fail,`
		830	`double abstol, double intol, void *ex, int type, int trace,`
7345	ripley	831	`int fncount, int grcount, int maxit)`
		832	`{`
10470	maechler	833	`Rboolean accpoint;`
7345	ripley	834	`double c, g, *t;`
		835	`int count, cycle, cyclimit;`
		836	`double f;`
		837	`double G1, G2, G3, gradproj;`
		838	`int funcount=0, gradcount=0, i;`
		839	`double newstep, oldstep, setstep, steplength=1.0;`
10229	ripley	840	`double tol;`
7345	ripley	841
18222	ripley	842	`if (maxit <= 0) {`
		843	`*Fmin = fminfn(n, Bvec, ex);`
		844	`fncount = grcount = 0;`
		845	`*fail = FALSE;`
		846	`return;`
		847	`}`
10196	maechler	848	`if (trace) {`
		849	`Rprintf(" Conjugate gradients function minimiser\n");`
7345	ripley	850	`switch (type) {`
10196	maechler	851	`case 1: Rprintf("Method: Fletcher Reeves\n"); break;`
		852	`case 2: Rprintf("Method: Polak Ribiere\n"); break;`
		853	`case 3: Rprintf("Method: Beale Sorenson\n"); break;`
7878	ripley	854	`default:`
		855	`error("unknown type in CG method of optim");`
7345	ripley	856	`}`
10196	maechler	857	`}`
		858	`c = vect(n); g = vect(n); t = vect(n);`
		859
		860	`setstep = 1.7;`
7345	ripley	861	`*fail = 0;`
		862	`cyclimit = n;`
		863	`tol = intol * n * sqrt(intol);`
		864
7510	ripley	865	`if (trace) Rprintf("tolerance used in gradient test=%g\n", tol);`
18016	ripley	866	`f = fminfn(n, Bvec, ex);`
7345	ripley	867	`if (!R_FINITE(f)) {`
		868	`error("Function cannot be evaluated at initial parameters");`
		869	`} else {`
		870	`*Fmin = f;`
		871	`funcount = 1;`
		872	`gradcount = 0;`
		873	`do {`
		874	`for (i = 0; i < n; i++) {`
		875	`t[i] = 0.0;`
		876	`c[i] = 0.0;`
		877	`}`
		878	`cycle = 0;`
		879	`oldstep = 1.0;`
		880	`count = 0;`
		881	`do {`
		882	`cycle++;`
7510	ripley	883	`if (trace) {`
7345	ripley	884	`Rprintf("%d %d %f\n", gradcount, funcount, *Fmin);`
		885	`Rprintf("parameters ");`
		886	`for (i = 1; i <= n; i++) {`
		887	`Rprintf("%10.5f ", Bvec[i - 1]);`
10196	maechler	888	`if (i / 7 * 7 == i && i < n)`
7345	ripley	889	`Rprintf("\n");`
		890	`}`
		891	`Rprintf("\n");`
		892	`}`
		893	`gradcount++;`
7510	ripley	894	`if (gradcount > maxit) {`
7345	ripley	895	`*fncount = funcount;`
		896	`*grcount = gradcount;`
		897	`*fail = 1;`
		898	`return;`
		899	`}`
18016	ripley	900	`fmingr(n, Bvec, g, ex);`
7345	ripley	901	`G1 = 0.0;`
		902	`G2 = 0.0;`
		903	`for (i = 0; i < n; i++) {`
		904	`X[i] = Bvec[i];`
		905	`switch (type) {`
		906
		907	`case 1: /* Fletcher-Reeves */`
10196	maechler	908	`G1 += g[i] * g[i];`
		909	`G2 += c[i] * c[i];`
7345	ripley	910	`break;`
		911
		912	`case 2: /* Polak-Ribiere */`
		913	`G1 += g[i] * (g[i] - c[i]);`
10196	maechler	914	`G2 += c[i] * c[i];`
7345	ripley	915	`break;`
		916
		917	`case 3: /* Beale-Sorenson */`
		918	`G1 += g[i] * (g[i] - c[i]);`
		919	`G2 += t[i] * (g[i] - c[i]);`
		920	`break;`
7878	ripley	921
		922	`default:`
		923	`error("unknown type in CG method of optim");`
7345	ripley	924	`}`
		925	`c[i] = g[i];`
		926	`}`
		927	`if (G1 > tol) {`
		928	`if (G2 > 0.0)`
		929	`G3 = G1 / G2;`
		930	`else`
		931	`G3 = 1.0;`
		932	`gradproj = 0.0;`
		933	`for (i = 0; i < n; i++) {`
		934	`t[i] = t[i] * G3 - g[i];`
		935	`gradproj += t[i] * g[i];`
		936	`}`
		937	`steplength = oldstep;`
		938
10470	maechler	939	`accpoint = FALSE;`
7345	ripley	940	`do {`
		941	`count = 0;`
		942	`for (i = 0; i < n; i++) {`
		943	`Bvec[i] = X[i] + steplength * t[i];`
		944	`if (reltest + X[i] == reltest + Bvec[i])`
		945	`count++;`
		946	`}`
		947	`if (count < n) {`
18016	ripley	948	`f = fminfn(n, Bvec, ex);`
7345	ripley	949	`funcount++;`
		950	`accpoint = (R_FINITE(f) &&`
		951	`f <= Fmin + gradproj steplength * acctol);`
		952
		953	`if (!accpoint) {`
		954	`steplength *= stepredn;`
7510	ripley	955	`if (trace) Rprintf("*");`
7345	ripley	956	`}`
		957	`}`
		958	`} while (!(count == n \|\| accpoint));`
		959	`if (count < n) {`
		960	`newstep = 2 * (f - Fmin - gradproj steplength);`
		961	`if (newstep > 0) {`
		962	`newstep = -(gradproj * steplength * steplength / newstep);`
		963	`for (i = 0; i < n; i++)`
		964	`Bvec[i] = X[i] + newstep * t[i];`
		965	`*Fmin = f;`
18016	ripley	966	`f = fminfn(n, Bvec, ex);`
7345	ripley	967	`funcount++;`
		968	`if (f < *Fmin) {`
		969	`*Fmin = f;`
7510	ripley	970	`if (trace) Rprintf(" i< ");`
7345	ripley	971	`} else {`
7510	ripley	972	`if (trace) Rprintf(" i> ");`
7345	ripley	973	`for (i = 0; i < n; i++)`
		974	`Bvec[i] = X[i] + steplength * t[i];`
		975	`}`
		976	`}`
		977	`}`
		978	`}`
		979	`oldstep = setstep * steplength;`
		980	`if (oldstep > 1.0)`
		981	`oldstep = 1.0;`
		982	`} while ((count != n) && (G1 > tol) && (cycle != cyclimit));`
		983
8186	ripley	984	`} while ((cycle != 1) \|\|`
7345	ripley	985	`((count != n) && (G1 > tol) && *Fmin > abstol));`
		986
		987	`}`
7510	ripley	988	`if (trace) {`
7345	ripley	989	`Rprintf("Exiting from conjugate gradients minimizer\n");`
		990	`Rprintf(" %d function evaluations used\n", funcount);`
		991	`Rprintf(" %d gradient evaluations used\n", gradcount);`
		992	`}`
		993	`*fncount = funcount;`
		994	`*grcount = gradcount;`
		995	`}`
7464	ripley	996
8186	ripley	997	`void lbfgsb(int n, int m, double x, double l, double u, int nbd,`
19233	maechler	998	`double Fmin, optimfn fminfn, optimgr fmingr, int fail,`
18016	ripley	999	`void *ex, double factr, double pgtol,`
10887	maechler	1000	`int fncount, int grcount, int maxit, char *msg,`
10853	ripley	1001	`int trace, int nREPORT)`
7464	ripley	1002	`{`
		1003	`char task[60];`
7505	ripley	1004	`double f, g, dsave[29], wa;`
10853	ripley	1005	`int tr = -1, iter = 0, *iwa, isave[44], lsave[4];`
7464	ripley	1006
10853	ripley	1007	`if (nREPORT <= 0)`
		1008	`error("REPORT must be > 0 (method = \"L-BFGS-B\")");`
		1009	`switch(trace) {`
		1010	`case 2: tr = 0; break;`
		1011	`case 3: tr = nREPORT; break;`
		1012	`case 4: tr = 99; break;`
		1013	`case 5: tr = 100; break;`
		1014	`case 6: tr = 101; break;`
		1015	`default: tr = -1; break;`
		1016	`}`
10887	maechler	1017
7464	ripley	1018	`*fail = 0;`
		1019	`g = vect(n);`
7505	ripley	1020	`wa = vect(2mn+4n+11mm+8m);`
		1021	`iwa = (int ) R_alloc(3n, sizeof(int));`
7464	ripley	1022	`strcpy(task, "START");`
		1023	`while(1) {`
13624	maechler	1024	`/* Main workhorse setulb() from ../appl/lbfgsb.c : */`
10887	maechler	1025	`setulb(n, m, x, l, u, nbd, &f, g, factr, &pgtol, wa, iwa, task,`
10853	ripley	1026	`tr, lsave, isave, dsave);`
7464	ripley	1027	`/* Rprintf("in lbfgsb - %s\n", task);*/`
		1028	`if (strncmp(task, "FG", 2) == 0) {`
18016	ripley	1029	`f = fminfn(n, x, ex);`
8186	ripley	1030	`if (!R_FINITE(f))`
7975	ripley	1031	`error("L-BFGS-B needs finite values of fn");`
18016	ripley	1032	`fmingr(n, x, g, ex);`
7464	ripley	1033	`} else if (strncmp(task, "NEW_X", 5) == 0) {`
10853	ripley	1034	`if(trace == 1 && (iter % nREPORT == 0)) {`
		1035	`Rprintf("iter %4d value %f\n", iter, f);`
		1036	`}`
7464	ripley	1037	`if (++iter > maxit) {`
		1038	`*fail = 1;`
		1039	`break;`
		1040	`}`
		1041	`} else if (strncmp(task, "WARN", 4) == 0) {`
		1042	`*fail = 51;`
		1043	`break;`
		1044	`} else if (strncmp(task, "CONV", 4) == 0) {`
		1045	`break;`
		1046	`} else if (strncmp(task, "ERROR", 5) == 0) {`
		1047	`*fail = 52;`
		1048	`break;`
10853	ripley	1049	`} else { /* some other condition that is not supposed to happen */`
		1050	`*fail = 52;`
10887	maechler	1051	`break;`
7464	ripley	1052	`}`
		1053	`}`
		1054	`*Fmin = f;`
		1055	`fncount = grcount = isave[33];`
10853	ripley	1056	`if (trace) {`
		1057	`Rprintf("final value %f \n", *Fmin);`
		1058	`if (iter < maxit && *fail == 0) Rprintf("converged\n");`
		1059	`else Rprintf("stopped after %i iterations\n", iter);`
10887	maechler	1060	`}`
7464	ripley	1061	`strcpy(msg, task);`
		1062	`}`
		1063
		1064
8186	ripley	1065	`#define E1 1.7182818 /* exp(1.0)-1.0 */`
		1066	`#define STEPS 100`
		1067
18031	ripley	1068	`void samin(int n, double pb, double yb, optimfn fminfn, int maxit,`
		1069	`int tmax, double ti, int trace, void *ex)`
8186	ripley	1070
		1071	`/* Given a starting point pb[0..n-1], simulated annealing minimization`
		1072	`is performed on the function fminfn. The starting temperature`
		1073	`is input as ti. To make sann work silently set trace to zero.`
		1074	`sann makes in total maxit function evaluations, tmax`
		1075	`evaluations at each temperature. Returned quantities are pb`
		1076	`(the location of the minimum), and yb (the minimum value of`
		1077	`the function func). Author: Adrian Trapletti`
		1078	`*/`
		1079	`{`
21689	ripley	1080	`long j;`
8186	ripley	1081	`int k, its, itdoc;`
		1082	`double t, y, dy, ytry, scale;`
		1083	`double p, dp, *ptry;`
		1084
		1085	`p = vect (n); dp = vect (n); ptry = vect (n);`
		1086	`GetRNGstate();`
18016	ripley	1087	`yb = fminfn (n, pb, ex); / init best system state pb, yb /`
8186	ripley	1088	`if (!R_FINITE(yb)) yb = big;`
		1089	`for (j = 0; j < n; j++) p[j] = pb[j];`
		1090	`y = yb; / init system state p, y */`
		1091	`if (trace)`
		1092	`{`
		1093	`Rprintf ("sann objective function values\n");`
		1094	`Rprintf ("initial value %f\n", *yb);`
		1095	`}`
		1096	`scale = 1.0/ti;`
		1097	`its = itdoc = 1;`
21689	ripley	1098	`while (its < maxit) { /* cool down system */`
8186	ripley	1099	`t = ti/log((double)its + E1); /* temperature annealing schedule */`
		1100	`k = 1;`
		1101	`while ((k <= tmax) && (its < maxit)) /* iterate at constant temperature */`
		1102	`{`
21689	ripley	1103	`genptry(n, p, ptry, scale * t, ex); /* generate new candidate point */`
18016	ripley	1104	`ytry = fminfn (n, ptry, ex);`
8186	ripley	1105	`if (!R_FINITE(ytry)) ytry = big;`
		1106	`dy = ytry - y;`
		1107	`if ((dy <= 0.0) \|\| (unif_rand() < exp(-dy/t))) { /* accept new point? */`
		1108	`for (j = 0; j < n; j++) p[j] = ptry[j];`
		1109	`y = ytry; /* update system state p, y */`
		1110	`if (y <= yb) / if system state is best, then update best system state pb, yb /`
		1111	`{`
		1112	`for (j = 0; j < n; j++) pb[j] = p[j];`
		1113	`*yb = y;`
		1114	`}`
		1115	`}`
		1116	`its++; k++;`
		1117	`}`
		1118	`if ((trace) && ((itdoc % STEPS) == 0))`
		1119	`Rprintf("iter %8d value %f\n", its - 1, *yb);`
		1120	`itdoc++;`
		1121	`}`
		1122	`if (trace)`
		1123	`{`
		1124	`Rprintf ("final value %f\n", *yb);`
		1125	`Rprintf ("sann stopped after %d iterations\n", its - 1);`
		1126	`}`
		1127	`PutRNGstate();`
		1128	`}`
		1129
		1130	`#undef E1`
		1131	`#undef STEPS`

The R Project SVN R

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