1 /*
2 * Routines used in calculating least squares solutions in a
3 * nonlinear model in nls library for R.
4 *
5 * Copyright 2005-2020 The R Core Team
6 * Copyright 2006 The R Foundation
7 * Copyright 1999-2001 Douglas M. Bates
8 * Saikat DebRoy
9 *
10 *
11 * This program is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation; either version 2 of the License, or
14 * (at your option) any later version.
15 *
16 * This program is distributed in the hope that it will be
17 * useful, but WITHOUT ANY WARRANTY; without even the implied
18 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
19 * PURPOSE. See the GNU General Public License for more
20 * details.
21 *
22 * You should have received a copy of the GNU General Public
23 * License along with this program; if not, a copy is available at
24 * https://www.R-project.org/Licenses/
25 */
26
27 #include <stdlib.h>
28 #include <string.h>
29 #include <math.h>
30 #include <float.h>
31 #include <R.h>
32 #include <Rinternals.h>
33 #include "nls.h"
34
35 #ifndef MIN
36 #define MIN(a,b) (((a)<(b))?(a):(b))
37 #endif
38
39 /*
40 * get the list element named str. names is the name attribute of list
41 */
42
43 static SEXP
getListElement(SEXP list,SEXP names,const char * str)44 getListElement(SEXP list, SEXP names, const char *str)
45 {
46 SEXP elmt = (SEXP) NULL;
47 const char *tempChar;
48 int i;
49
50 for (i = 0; i < length(list); i++) {
51 tempChar = CHAR(STRING_ELT(names, i)); /* ASCII only */
52 if( strcmp(tempChar,str) == 0) {
53 elmt = VECTOR_ELT(list, i);
54 break;
55 }
56 }
57 return elmt;
58 }
59
60 /*
61 * put some convergence-related information into list
62 */
63 static SEXP
ConvInfoMsg(char * msg,int iter,int whystop,double fac,double minFac,int maxIter,double convNew)64 ConvInfoMsg(char* msg, int iter, int whystop, double fac,
65 double minFac, int maxIter, double convNew)
66 {
67 const char *nms[] = {"isConv", "finIter", "finTol",
68 "stopCode", "stopMessage", ""};
69 SEXP ans;
70 PROTECT(ans = mkNamed(VECSXP, nms));
71
72 SET_VECTOR_ELT(ans, 0, ScalarLogical(whystop == 0)); /* isConv */
73 SET_VECTOR_ELT(ans, 1, ScalarInteger(iter)); /* finIter */
74 SET_VECTOR_ELT(ans, 2, ScalarReal (convNew)); /* finTol */
75 SET_VECTOR_ELT(ans, 3, ScalarInteger(whystop)); /* stopCode */
76 SET_VECTOR_ELT(ans, 4, mkString(msg)); /* stopMessage */
77
78 UNPROTECT(1);
79 return ans;
80 }
81
82
83 /*
84 * call to nls_iter from R --- .Call("nls_iter", m, control, doTrace)
85 * where m and control are nlsModel and nlsControl objects
86 * doTrace is a logical value.
87 * m is modified; the return value is a "convergence-information" list.
88 */
89 SEXP
nls_iter(SEXP m,SEXP control,SEXP doTraceArg)90 nls_iter(SEXP m, SEXP control, SEXP doTraceArg)
91 {
92 int doTrace = asLogical(doTraceArg);
93
94 if(!isNewList(control))
95 error(_("'control' must be a list"));
96 if(!isNewList(m))
97 error(_("'m' must be a list"));
98
99 SEXP tmp, conv;
100 PROTECT(tmp = getAttrib(control, R_NamesSymbol));
101
102 conv = getListElement(control, tmp, "maxiter");
103 if(conv == NULL || !isNumeric(conv))
104 error(_("'%s' absent"), "control$maxiter");
105 int maxIter = asInteger(conv);
106
107 conv = getListElement(control, tmp, "tol");
108 if(conv == NULL || !isNumeric(conv))
109 error(_("'%s' absent"), "control$tol");
110 double tolerance = asReal(conv);
111
112 conv = getListElement(control, tmp, "minFactor");
113 if(conv == NULL || !isNumeric(conv))
114 error(_("'%s' absent"), "control$minFactor");
115 double minFac = asReal(conv);
116
117 conv = getListElement(control, tmp, "warnOnly");
118 if(conv == NULL || !isLogical(conv))
119 error(_("'%s' absent"), "control$warnOnly");
120 int warnOnly = asLogical(conv);
121
122 conv = getListElement(control, tmp, "printEval");
123 if(conv == NULL || !isLogical(conv))
124 error(_("'%s' absent"), "control$printEval");
125 Rboolean printEval = asLogical(conv);
126
127 // now get parts from 'm' ---------------------------------
128 tmp = getAttrib(m, R_NamesSymbol);
129
130 conv = getListElement(m, tmp, "conv");
131 if(conv == NULL || !isFunction(conv))
132 error(_("'%s' absent"), "m$conv()");
133 PROTECT(conv = lang1(conv));
134
135 SEXP incr = getListElement(m, tmp, "incr");
136 if(incr == NULL || !isFunction(incr))
137 error(_("'%s' absent"), "m$incr()");
138 PROTECT(incr = lang1(incr));
139
140 SEXP deviance = getListElement(m, tmp, "deviance");
141 if(deviance == NULL || !isFunction(deviance))
142 error(_("'%s' absent"), "m$deviance()");
143 PROTECT(deviance = lang1(deviance));
144
145 SEXP trace = getListElement(m, tmp, "trace");
146 if(trace == NULL || !isFunction(trace))
147 error(_("'%s' absent"), "m$trace()");
148 PROTECT(trace = lang1(trace));
149
150 SEXP setPars = getListElement(m, tmp, "setPars");
151 if(setPars == NULL || !isFunction(setPars))
152 error(_("'%s' absent"), "m$setPars()");
153 PROTECT(setPars);
154
155 SEXP getPars = getListElement(m, tmp, "getPars");
156 if(getPars == NULL || !isFunction(getPars))
157 error(_("'%s' absent"), "m$getPars()");
158 PROTECT(getPars = lang1(getPars));
159
160 SEXP pars = PROTECT(eval(getPars, R_GlobalEnv));
161 int nPars = LENGTH(pars);
162
163 double dev = asReal(eval(deviance, R_GlobalEnv));
164 if(doTrace) eval(trace,R_GlobalEnv);
165
166 double fac = 1.0;
167 Rboolean hasConverged = FALSE;
168 SEXP newPars = PROTECT(allocVector(REALSXP, nPars));
169 int evaltotCnt = 1;
170 double convNew = -1. /* -Wall */;
171 int i;
172 #define CONV_INFO_MSG(_STR_, _I_) \
173 ConvInfoMsg(_STR_, i, _I_, fac, minFac, maxIter, convNew)
174
175 #define NON_CONV_FINIS(_ID_, _MSG_) \
176 if(warnOnly) { \
177 warning(_MSG_); \
178 return CONV_INFO_MSG(_MSG_, _ID_); \
179 } \
180 else \
181 error(_MSG_);
182
183 #define NON_CONV_FINIS_1(_ID_, _MSG_, _A1_) \
184 if(warnOnly) { \
185 char msgbuf[1000]; \
186 warning(_MSG_, _A1_); \
187 snprintf(msgbuf, 1000, _MSG_, _A1_); \
188 return CONV_INFO_MSG(msgbuf, _ID_); \
189 } \
190 else \
191 error(_MSG_, _A1_);
192
193 #define NON_CONV_FINIS_2(_ID_, _MSG_, _A1_, _A2_) \
194 if(warnOnly) { \
195 char msgbuf[1000]; \
196 warning(_MSG_, _A1_, _A2_); \
197 snprintf(msgbuf, 1000, _MSG_, _A1_, _A2_); \
198 return CONV_INFO_MSG(msgbuf, _ID_); \
199 } \
200 else \
201 error(_MSG_, _A1_, _A2_);
202
203 for (i = 0; i < maxIter; i++) { // ---------------------------------------------
204
205 if((convNew = asReal(eval(conv, R_GlobalEnv))) <= tolerance) {
206 hasConverged = TRUE;
207 break;
208 }
209
210 SEXP newIncr = PROTECT(eval(incr, R_GlobalEnv));
211 double
212 *par = REAL(pars),
213 *npar = REAL(newPars),
214 *nIncr = REAL(newIncr);
215 int evalCnt = -1;
216 if(printEval)
217 evalCnt = 1;
218
219 while(fac >= minFac) { // 1-dim "line search"
220 if(printEval) {
221 Rprintf(" It. %3d, fac= %11.6g, eval (no.,total): (%2d,%3d):",
222 i+1, fac, evalCnt, evaltotCnt);
223 evalCnt++;
224 evaltotCnt++;
225 }
226 for(int j = 0; j < nPars; j++)
227 npar[j] = par[j] + fac * nIncr[j];
228
229 PROTECT(tmp = lang2(setPars, newPars));
230 if (asLogical(eval(tmp, R_GlobalEnv))) { /* singular gradient */
231 UNPROTECT(11);
232
233 NON_CONV_FINIS(1, _("singular gradient"));
234 }
235 UNPROTECT(1);
236
237 double newDev = asReal(eval(deviance, R_GlobalEnv));
238 if(printEval)
239 Rprintf(" new dev = %g\n", newDev);
240 if(newDev <= dev) {
241 dev = newDev;
242 tmp = newPars;
243 newPars = pars;
244 pars = tmp;
245 fac = MIN(2*fac, 1);
246 break;
247 } // else
248 fac /= 2.;
249 }
250 UNPROTECT(1);
251 if(doTrace) eval(trace, R_GlobalEnv);
252 if( fac < minFac ) {
253 UNPROTECT(9);
254 NON_CONV_FINIS_2(2,
255 _("step factor %g reduced below 'minFactor' of %g"),
256 fac, minFac);
257 }
258 }
259
260 UNPROTECT(9);
261 if(!hasConverged) {
262 NON_CONV_FINIS_1(3,
263 _("number of iterations exceeded maximum of %d"),
264 maxIter);
265 }
266 else
267 return CONV_INFO_MSG(_("converged"), 0);
268 }
269 #undef CONV_INFO_MSG
270 #undef NON_CONV_FINIS
271 #undef NON_CONV_FINIS_1
272 #undef NON_CONV_FINIS_2
273
274
275 /*
276 * call to numeric_deriv from R -
277 * .Call("numeric_deriv", expr, theta, rho, dir = 1., eps = .Machine$double.eps, central=FALSE)
278 * Returns: ans
279 */
280 SEXP
numeric_deriv(SEXP expr,SEXP theta,SEXP rho,SEXP dir,SEXP eps_,SEXP centr)281 numeric_deriv(SEXP expr, SEXP theta, SEXP rho, SEXP dir, SEXP eps_, SEXP centr)
282 {
283 if(!isString(theta))
284 error(_("'theta' should be of type character"));
285 if (isNull(rho)) {
286 error(_("use of NULL environment is defunct"));
287 rho = R_BaseEnv;
288 } else
289 if(!isEnvironment(rho))
290 error(_("'rho' should be an environment"));
291 int nprot = 3;
292 if(TYPEOF(dir) != REALSXP) {
293 PROTECT(dir = coerceVector(dir, REALSXP)); nprot++;
294 }
295 if(LENGTH(dir) != LENGTH(theta))
296 error(_("'dir' is not a numeric vector of the correct length"));
297 Rboolean central = asLogical(centr);
298 if(central == NA_LOGICAL)
299 error(_("'central' is NA, but must be TRUE or FALSE"));
300 SEXP rho1 = PROTECT(R_NewEnv(rho, FALSE, 0));
301 nprot++;
302 SEXP
303 pars = PROTECT(allocVector(VECSXP, LENGTH(theta))),
304 ans = PROTECT(duplicate(eval(expr, rho1)));
305 double *rDir = REAL(dir), *res = NULL; // -Wall
306 #define CHECK_FN_VAL(_r_, _ANS_) do { \
307 if(!isReal(_ANS_)) { \
308 SEXP temp = coerceVector(_ANS_, REALSXP); \
309 UNPROTECT(1);/*: _ANS_ *must* have been the last PROTECT() ! */ \
310 PROTECT(_ANS_ = temp); \
311 } \
312 _r_ = REAL(_ANS_); \
313 for(int i = 0; i < LENGTH(_ANS_); i++) { \
314 if (!R_FINITE(_r_[i])) \
315 error(_("Missing value or an infinity produced when evaluating the model")); \
316 } \
317 } while(0)
318
319 CHECK_FN_VAL(res, ans);
320
321 const void *vmax = vmaxget();
322 int lengthTheta = 0;
323 for(int i = 0; i < LENGTH(theta); i++) {
324 const char *name = translateChar(STRING_ELT(theta, i));
325 SEXP s_name = install(name);
326 SEXP temp = findVar(s_name, rho1);
327 if(isInteger(temp))
328 error(_("variable '%s' is integer, not numeric"), name);
329 if(!isReal(temp))
330 error(_("variable '%s' is not numeric"), name);
331 // We'll be modifying the variable, so need to make a copy PR#15849
332 defineVar(s_name, temp = duplicate(temp), rho1);
333 MARK_NOT_MUTABLE(temp);
334 SET_VECTOR_ELT(pars, i, temp);
335 lengthTheta += LENGTH(VECTOR_ELT(pars, i));
336 }
337 vmaxset(vmax);
338 SEXP gradient = PROTECT(allocMatrix(REALSXP, LENGTH(ans), lengthTheta));
339 double *grad = REAL(gradient);
340 double eps = asReal(eps_); // was hardcoded sqrt(DOUBLE_EPS) { ~= 1.49e-08, typically}
341 for(int start = 0, i = 0; i < LENGTH(theta); i++) {
342 double *pars_i = REAL(VECTOR_ELT(pars, i));
343 for(int j = 0; j < LENGTH(VECTOR_ELT(pars, i)); j++, start += LENGTH(ans)) {
344 double
345 origPar = pars_i[j],
346 xx = fabs(origPar),
347 delta = (xx == 0) ? eps : xx*eps;
348 pars_i[j] += rDir[i] * delta;
349 SEXP ans_del = PROTECT(eval(expr, rho1));
350 double *rDel = NULL;
351 CHECK_FN_VAL(rDel, ans_del);
352 if(central) {
353 pars_i[j] = origPar - rDir[i] * delta;
354 SEXP ans_de2 = PROTECT(eval(expr, rho1));
355 double *rD2 = NULL;
356 CHECK_FN_VAL(rD2, ans_de2);
357 for(int k = 0; k < LENGTH(ans); k++) {
358 grad[start + k] = rDir[i] * (rDel[k] - rD2[k])/(2 * delta);
359 }
360 } else { // forward difference (previously hardwired):
361 for(int k = 0; k < LENGTH(ans); k++) {
362 grad[start + k] = rDir[i] * (rDel[k] - res[k])/delta;
363 }
364 }
365 UNPROTECT(central ? 2 : 1); // ansDel & possibly ans
366 pars_i[j] = origPar;
367 }
368 }
369 setAttrib(ans, install("gradient"), gradient);
370 UNPROTECT(nprot);
371 return ans;
372 }
373