1 /*
2 * Copyright (C) 1998--2020 The R Core Team
3 *
4 * The authors of this software are Cleveland, Grosse, and Shyu.
5 * Copyright (c) 1989, 1992 by AT&T.
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose without fee is hereby granted, provided that this entire notice
8 * is included in all copies of any software which is or includes a copy
9 * or modification of this software and in all copies of the supporting
10 * documentation for such software.
11 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
12 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
13 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
14 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
15 */
16
17 /* <UTF8> chars are handled as whole strings.
18 They are passed from Fortran so had better be ASCII.
19 */
20
21 /*
22 * Altered by B.D. Ripley to use F77_*, declare routines before use.
23 *
24 * 'protoize'd to ANSI C headers; indented: M.Maechler
25 *
26 * Changes to the C/Fortran interface to support LTO in May 2019.
27 *
28 * lowesd() : eliminate 'version' argument { also in ./loessf.f }
29 */
30
31 #ifdef HAVE_CONFIG_H
32 // for FC_LEN_T
33 #include <config.h>
34 #endif
35
36 #include <string.h>
37 #include <stdio.h>
38 #include <math.h>
39 #include <limits.h>
40
41 #include <Rmath.h> // R_pow_di()
42 #include "modreg.h"
43
44
45 /* Forward declarations */
46 static
47 void loess_workspace(int D, int N, double span, int degree,
48 int nonparametric, const int drop_square[],
49 int sum_drop_sqr, Rboolean setLf);
50
51 static
52 void loess_prune(int *parameter, int *a,
53 double *xi, double *vert, double *vval);
54 static
55 void loess_grow (int *parameter, int *a,
56 double *xi, double *vert, double *vval);
57
58 /* These (and many more) are in ./loessf.f : */
59 void F77_NAME(lowesa)(double*, int*, int*, int*, int*, double*, double*);
60 void F77_NAME(lowesb)(double*, double*, double*, double*, int*, int*, double*);
61 void F77_NAME(lowesc)(int*, double*, double*, double*, double*, double*);
62 void F77_NAME(lowesd)(int* iv, int* liv, int* lv, double *v,
63 int* d, int* n, double* f,
64 int* ideg, int* nf, int* nvmax, int* setlf);
65 void F77_NAME(lowese)(int*, double*, int*, double*, double*);
66 void F77_NAME(lowesf)(double*, double*, double*, int*, double*,
67 int*, double*, double*, int*, double*);
68 void F77_NAME(lowesl)(int*, double*, int*, double*, double*);
69 void F77_NAME(ehg169)(int*, int*, int*, int*, int*, int*,
70 double*, int*, double*, int*, int*, int*);
71 void F77_NAME(ehg196)(int*, int*, double*, double*);
72 /* exported (for loessf.f) : */
73 void F77_SUB(loesswarn)(int *i);
74 #ifdef FC_LEN_T
75 # include <stddef.h>
76 void F77_SUB(ehg183a)(char *s, int *nc,int *i,int *n,int *inc, FC_LEN_T c1);
77 void F77_SUB(ehg184a)(char *s, int *nc, double *x, int *n, int *inc, FC_LEN_T c1);
78 #else
79 void F77_SUB(ehg183a)(char *s, int *nc,int *i,int *n,int *inc);
80 void F77_SUB(ehg184a)(char *s, int *nc, double *x, int *n, int *inc);
81 #endif
82
83
84 #undef min
85 #undef max
86
87 #define min(x,y) ((x) < (y) ? (x) : (y))
88 #define max(x,y) ((x) > (y) ? (x) : (y))
89 #define GAUSSIAN 1
90 #define SYMMETRIC 0
91
92 // Global variables :
93 static int *iv = NULL, liv, lv, tau;
94 static double *v = NULL;
95
96 /* these are set in an earlier call to loess_workspace or loess_grow */
loess_free(void)97 static void loess_free(void)
98 {
99 Free(v);
100 Free(iv);
101 }
102
103 void
loess_raw(double * y,double * x,double * weights,double * robust,int * d,int * n,double * span,int * degree,int * nonparametric,int * drop_square,int * sum_drop_sqr,double * cell,char ** surf_stat,double * surface,int * parameter,int * a,double * xi,double * vert,double * vval,double * diagonal,double * trL,double * one_delta,double * two_delta,int * setLf)104 loess_raw(double *y, double *x, double *weights, double *robust, int *d,
105 int *n, double *span, int *degree, int *nonparametric,
106 int *drop_square, int *sum_drop_sqr, double *cell,
107 char **surf_stat, double *surface, int *parameter,
108 int *a, double *xi, double *vert, double *vval, double *diagonal,
109 double *trL, double *one_delta, double *two_delta, int *setLf)
110 {
111 int i0 = 0, one = 1, two = 2, nsing, i, k;
112 double *hat_matrix, *LL, d0=0.0;
113
114 *trL = 0;
115
116 loess_workspace(*d, *n, *span, *degree, *nonparametric, drop_square, *sum_drop_sqr, *setLf);
117 v[1] = *cell;/* = v(2) in Fortran (!) */
118
119 /* NB: surf_stat = (surface / statistics);
120 * statistics = "none" for all robustness iterations
121 */
122 if(!strcmp(*surf_stat, "interpolate/none")) { // default for loess.smooth() and robustness iter.
123 F77_CALL(lowesb)(x, y, robust, &d0, &i0, iv, v);
124 F77_CALL(lowese)(iv, v, n, x, surface);
125 loess_prune(parameter, a, xi, vert, vval);
126 }
127 else if (!strcmp(*surf_stat, "direct/none")) {
128 F77_CALL(lowesf)(x, y, robust, iv, v, n, x, &d0, &i0, surface);
129 }
130 else if (!strcmp(*surf_stat, "interpolate/1.approx")) { // default (trace.hat is "exact")
131 F77_CALL(lowesb)(x, y, weights, diagonal, &one, iv, v);
132 F77_CALL(lowese)(iv, v, n, x, surface);
133 nsing = iv[29];
134 for(i = 0; i < (*n); i++) *trL = *trL + diagonal[i];
135 F77_CALL(lowesa)(trL, n, d, &tau, &nsing, one_delta, two_delta);
136 loess_prune(parameter, a, xi, vert, vval);
137 }
138 else if (!strcmp(*surf_stat, "interpolate/2.approx")) { // default for trace.hat = "approximate"
139 // vvvvvvv (had 'robust' in R <= 3.2.x)
140 F77_CALL(lowesb)(x, y, weights, &d0, &i0, iv, v);
141 F77_CALL(lowese)(iv, v, n, x, surface);
142 nsing = iv[29];
143 F77_CALL(ehg196)(&tau, d, span, trL);
144 F77_CALL(lowesa)(trL, n, d, &tau, &nsing, one_delta, two_delta);
145 loess_prune(parameter, a, xi, vert, vval);
146 }
147 else if (!strcmp(*surf_stat, "direct/approximate")) {
148 F77_CALL(lowesf)(x, y, weights, iv, v, n, x, diagonal, &one, surface);
149 nsing = iv[29];
150 for(i = 0; i < (*n); i++) *trL = *trL + diagonal[i];
151 F77_CALL(lowesa)(trL, n, d, &tau, &nsing, one_delta, two_delta);
152 }
153 else if (!strcmp(*surf_stat, "interpolate/exact")) {
154 hat_matrix = (double *) R_alloc((*n)*(*n), sizeof(double));
155 LL = (double *) R_alloc((*n)*(*n), sizeof(double));
156 F77_CALL(lowesb)(x, y, weights, diagonal, &one, iv, v);
157 F77_CALL(lowesl)(iv, v, n, x, hat_matrix);
158 F77_CALL(lowesc)(n, hat_matrix, LL, trL, one_delta, two_delta);
159 F77_CALL(lowese)(iv, v, n, x, surface);
160 loess_prune(parameter, a, xi, vert, vval);
161 }
162 else if (!strcmp(*surf_stat, "direct/exact")) {
163 hat_matrix = (double *) R_alloc((*n)*(*n), sizeof(double));
164 LL = (double *) R_alloc((*n)*(*n), sizeof(double));
165 F77_CALL(lowesf)(x, y, weights, iv, v, n, x, hat_matrix, &two, surface);
166 F77_CALL(lowesc)(n, hat_matrix, LL, trL, one_delta, two_delta);
167 k = (*n) + 1;
168 for(i = 0; i < (*n); i++)
169 diagonal[i] = hat_matrix[i * k];
170 }
171 loess_free();
172 }
173
174 void
loess_dfit(double * y,double * x,double * x_evaluate,double * weights,double * span,int * degree,int * nonparametric,int * drop_square,int * sum_drop_sqr,int * d,int * n,int * m,double * fit)175 loess_dfit(double *y, double *x, double *x_evaluate, double *weights,
176 double *span, int *degree, int *nonparametric,
177 int *drop_square, int *sum_drop_sqr,
178 int *d, int *n, int *m, double *fit)
179 {
180 int i0 = 0;
181 double d0 = 0.0;
182
183 loess_workspace(*d, *n, *span, *degree, *nonparametric, drop_square, *sum_drop_sqr, FALSE);
184 F77_CALL(lowesf)(x, y, weights, iv, v, m, x_evaluate, &d0, &i0, fit);
185 loess_free();
186 }
187
188 void
loess_dfitse(double * y,double * x,double * x_evaluate,double * weights,double * robust,int * family,double * span,int * degree,int * nonparametric,int * drop_square,int * sum_drop_sqr,int * d,int * n,int * m,double * fit,double * L)189 loess_dfitse(double *y, double *x, double *x_evaluate, double *weights,
190 double *robust, int *family, double *span, int *degree,
191 int *nonparametric, int *drop_square,
192 int *sum_drop_sqr,
193 int *d, int *n, int *m, double *fit, double *L)
194 {
195 loess_workspace(*d, *n, *span, *degree, *nonparametric, drop_square, *sum_drop_sqr, FALSE);
196
197 int i2 = 2;
198 if(*family == GAUSSIAN)
199 F77_CALL(lowesf)(x, y, weights, iv, v, m, x_evaluate, L, &i2, fit);
200 else if(*family == SYMMETRIC)
201 {
202 int i0 = 0; double d0 = 0.0;
203 F77_CALL(lowesf)(x, y, weights, iv, v, m, x_evaluate, L, &i2, fit);
204 F77_CALL(lowesf)(x, y, robust, iv, v, m, x_evaluate, &d0, &i0, fit);
205 }
206 loess_free();
207 }
208
209 void
loess_ifit(int * parameter,int * a,double * xi,double * vert,double * vval,int * m,double * x_evaluate,double * fit)210 loess_ifit(int *parameter, int *a, double *xi, double *vert,
211 double *vval, int *m, double *x_evaluate, double *fit)
212 {
213 loess_grow(parameter, a, xi, vert, vval);
214 F77_CALL(lowese)(iv, v, m, x_evaluate, fit);
215 loess_free();
216 }
217
218 // Called from R's predLoess() when 'se = TRUE' (and the default surface == "interpolate")
219 void
loess_ise(double * y,double * x,double * x_evaluate,double * weights,double * span,int * degree,int * nonparametric,int * drop_square,int * sum_drop_sqr,double * cell,int * d,int * n,int * m,double * fit,double * L)220 loess_ise(double *y, double *x, double *x_evaluate, double *weights,
221 double *span, int *degree, int *nonparametric,
222 int *drop_square, int *sum_drop_sqr, double *cell,
223 int *d, int *n, int *m, double *fit, double *L)
224 {
225 loess_workspace(*d, *n, *span, *degree, *nonparametric, drop_square, *sum_drop_sqr, TRUE);
226
227 int i0 = 0; double d0 = 0.0;
228 v[1] = *cell;
229 F77_CALL(lowesb)(x, y, weights, &d0, &i0, iv, v);
230 F77_CALL(lowesl)(iv, v, m, x_evaluate, L);
231 loess_free();
232 }
233
234 // Set global variables tau, lv, liv , and allocate global arrays v[1..lv], iv[1..liv]
235 void
loess_workspace(int D,int N,double span,int degree,int nonparametric,const int drop_square[],int sum_drop_sqr,Rboolean setLf)236 loess_workspace(int D, int N, double span, int degree,
237 int nonparametric, const int drop_square[],
238 int sum_drop_sqr, Rboolean setLf)
239 {
240 int nvmax = max(200, N),
241 nf = min(N, (int) floor(N * span + 1e-5));
242 if(nf <= 0) error(_("span is too small"));
243 // NB: D := ncol(x) is <= 3
244 int tau0 = (degree > 1) ? ((D + 2) * (D + 1)) / 2 : (D + 1);
245 tau = tau0 - sum_drop_sqr;
246 double dlv = 50 + (3 * D + 3) * (double)nvmax + N + (tau0 + 2.) * nf;
247 double dliv = 50 + (R_pow_di(2., D) + 4.) * nvmax + 2. * N;
248 if(setLf) {
249 dlv += (D + 1.) * nf * (double)nvmax;
250 dliv += nf * (double)nvmax;
251 }
252
253 if (dlv < INT_MAX && dliv < INT_MAX) { // set the global vars
254 lv = (int) dlv;
255 liv = (int) dliv;
256 } else {
257 error(_("workspace required (%.0f) is too large%s."), max(dlv, dliv),
258 setLf ? _(" probably because of setting 'se = TRUE'") : "");
259 }
260
261 iv = Calloc(liv, int);
262 v = Calloc(lv, double);
263
264 F77_CALL(lowesd)(iv, &liv, &lv, v, &D, &N, &span,
265 °ree, &nf, &nvmax, (int *) &setLf);
266 iv[32] = nonparametric;
267 for(int i = 0; i < D; i++)
268 iv[i + 40] = drop_square[i];
269 }
270
271 static void
loess_prune(int * parameter,int * a,double * xi,double * vert,double * vval)272 loess_prune(int *parameter, int *a, double *xi, double *vert,
273 double *vval)
274 {
275 int d, vc, a1, v1, xi1, vv1, nc, nv, nvmax, i, k;
276
277 d = iv[1];
278 vc = iv[3] - 1;
279 nc = iv[4];
280 nv = iv[5];
281 a1 = iv[6] - 1;
282 v1 = iv[10] - 1;
283 xi1 = iv[11] - 1;
284 vv1 = iv[12] - 1;
285 nvmax = iv[13];
286
287 for(i = 0; i < 5; i++)
288 parameter[i] = iv[i + 1];
289 parameter[5] = iv[21] - 1;
290 parameter[6] = iv[14] - 1;
291
292 for(i = 0; i < d; i++) {
293 k = nvmax * i;
294 vert[i] = v[v1 + k];
295 vert[i + d] = v[v1 + vc + k];
296 }
297 for(i = 0; i < nc; i++) {
298 xi[i] = v[xi1 + i];
299 a[i] = iv[a1 + i];
300 }
301 k = (d + 1) * nv;
302 for(i = 0; i < k; i++)
303 vval[i] = v[vv1 + i];
304 }
305
306 static void
loess_grow(int * parameter,int * a,double * xi,double * vert,double * vval)307 loess_grow(int *parameter, int *a, double *xi,
308 double *vert, double *vval)
309 {
310 int d, vc, nc, nv, a1, v1, xi1, vv1, i, k;
311
312 d = parameter[0];
313 vc = parameter[2];
314 nc = parameter[3];
315 nv = parameter[4];
316 liv = parameter[5];
317 lv = parameter[6];
318 iv = Calloc(liv, int);
319 v = Calloc(lv, double);
320
321 iv[1] = d;
322 iv[2] = parameter[1];
323 iv[3] = vc;
324 iv[5] = iv[13] = nv;
325 iv[4] = iv[16] = nc;
326 iv[6] = 50;
327 iv[7] = iv[6] + nc;
328 iv[8] = iv[7] + vc * nc;
329 iv[9] = iv[8] + nc;
330 iv[10] = 50;
331 iv[12] = iv[10] + nv * d;
332 iv[11] = iv[12] + (d + 1) * nv;
333 iv[27] = 173;
334
335 v1 = iv[10] - 1;
336 xi1 = iv[11] - 1;
337 a1 = iv[6] - 1;
338 vv1 = iv[12] - 1;
339
340 for(i = 0; i < d; i++) {
341 k = nv * i;
342 v[v1 + k] = vert[i];
343 v[v1 + vc - 1 + k] = vert[i + d];
344 }
345 for(i = 0; i < nc; i++) {
346 v[xi1 + i] = xi[i];
347 iv[a1 + i] = a[i];
348 }
349 k = (d + 1) * nv;
350 for(i = 0; i < k; i++)
351 v[vv1 + i] = vval[i];
352
353 F77_CALL(ehg169)(&d, &vc, &nc, &nc, &nv, &nv, v+v1, iv+a1,
354 v+xi1, iv+iv[7]-1, iv+iv[8]-1, iv+iv[9]-1);
355 }
356
357
358 /* begin ehg's FORTRAN-callable C-codes */
359 #define MSG(_m_) msg = _(_m_) ; break ;
360
F77_SUB(loesswarn)361 void F77_SUB(loesswarn)(int *i)
362 {
363 char *msg, msg2[50];
364
365 switch(*i){
366 case 100:MSG("wrong version number in lowesd. Probably typo in caller.")
367 case 101:MSG("d>dMAX in ehg131. Need to recompile with increased dimensions.")
368 case 102:MSG("liv too small. (Discovered by lowesd)")
369 case 103:MSG("lv too small. (Discovered by lowesd)")
370 case 104:MSG("span too small. fewer data values than degrees of freedom.")
371 case 105:MSG("k>d2MAX in ehg136. Need to recompile with increased dimensions.")
372 case 106:MSG("lwork too small")
373 case 107:MSG("invalid value for kernel")
374 case 108:MSG("invalid value for ideg")
375 case 109:MSG("lowstt only applies when kernel=1.")
376 case 110:MSG("not enough extra workspace for robustness calculation")
377 case 120:MSG("zero-width neighborhood. make span bigger")
378 case 121:MSG("all data on boundary of neighborhood. make span bigger")
379 case 122:MSG("extrapolation not allowed with blending")
380 case 123:MSG("ihat=1 (diag L) in l2fit only makes sense if z=x (eval=data).")
381 case 171:MSG("lowesd must be called first.")
382 case 172:MSG("lowesf must not come between lowesb and lowese, lowesr, or lowesl.")
383 case 173:MSG("lowesb must come before lowese, lowesr, or lowesl.")
384 case 174:MSG("lowesb need not be called twice.")
385 case 175:MSG("need setLf=.true. for lowesl.")
386 case 180:MSG("nv>nvmax in cpvert.")
387 case 181:MSG("nt>20 in eval.")
388 case 182:MSG("svddc failed in l2fit.")
389 case 183:MSG("didnt find edge in vleaf.")
390 case 184:MSG("zero-width cell found in vleaf.")
391 case 185:MSG("trouble descending to leaf in vleaf.")
392 case 186:MSG("insufficient workspace for lowesf.")
393 case 187:MSG("insufficient stack space")
394 case 188:MSG("lv too small for computing explicit L")
395 case 191:MSG("computed trace L was negative; something is wrong!")
396 case 192:MSG("computed delta was negative; something is wrong!")
397 case 193:MSG("workspace in loread appears to be corrupted")
398 case 194:MSG("trouble in l2fit/l2tr")
399 case 195:MSG("only constant, linear, or quadratic local models allowed")
400 case 196:MSG("degree must be at least 1 for vertex influence matrix")
401 case 999:MSG("not yet implemented")
402 default: {
403 snprintf(msg2, 50, "Assert failed; error code %d\n",*i);
404 msg = msg2;
405 }
406 }
407 warning(msg);
408 }
409 #undef MSG
410
411 #ifdef FC_LEN_T
F77_SUB(ehg183a)412 void F77_SUB(ehg183a)(char *s, int *nc,int *i,int *n,int *inc, FC_LEN_T c1)
413 #else
414 void F77_SUB(ehg183a)(char *s, int *nc,int *i,int *n,int *inc)
415 #endif
416 {
417 int nnc = *nc;
418 char mess[4000], num[20];
419 strncpy(mess, s, nnc);
420 mess[nnc] = '\0';
421 for (int j = 0; j < *n; j++) {
422 snprintf(num, 20, " %d", i[j * *inc]);
423 strcat(mess, num);
424 }
425 strcat(mess,"\n");
426 warning(mess);
427 }
428
429 #ifdef FC_LEN_T
F77_SUB(ehg184a)430 void F77_SUB(ehg184a)(char *s, int *nc, double *x, int *n, int *inc, FC_LEN_T c1)
431 #else
432 void F77_SUB(ehg184a)(char *s, int *nc, double *x, int *n, int *inc)
433 #endif
434 {
435 int nnc = *nc;
436 char mess[4000], num[30];
437 strncpy(mess, s, nnc);
438 mess[nnc] = '\0';
439 for (int j = 0; j < *n; j++) {
440 snprintf(num, 30, " %.5g", x[j * *inc]);
441 strcat(mess, num);
442 }
443 strcat(mess,"\n");
444 warning(mess);
445 }
446