Pristine Start using Luke's original CLS 1.0 alpha 1
[tsl.git] / lib / lowess.c
blobfb3875fa6a61b00fbba48daf03cbc6ad6e02568a
1 /*
2 Translated from RATFOR lowess code of W. S. Cleveland as obtained from NETLIB
3 */
5 #include "xmath.h"
7 #define FALSE 0
8 #define TRUE 1
10 static double pow2(x) double x; { return(x * x); }
11 static double pow3(x) double x; { return(x * x * x); }
12 static double fmax(x,y) double x, y; { return (x > y ? x : y); }
14 lowess(x, y, n, f, nsteps, delta, ys, rw, res)
15 double *x, *y, f, delta, *ys, *rw, *res;
16 int n, nsteps;
18 int iter, ns, ok, nleft, nright, i, j, last, m1, m2;
19 double d1, d2, denom, alpha, cut, cmad, c9, c1, r;
21 if (n < 2) { ys[0] = y[0]; return(1); }
22 ns = max(min((int) (f * n), n), 2); /* at least two, at most n points */
23 for(iter = 1; iter <= nsteps + 1; iter++){ /* robustness iterations */
24 nleft = 0; nright = ns - 1;
25 last = -1; /* index of prev estimated point */
26 i = 0; /* index of current point */
27 do {
28 while(nright < n - 1){
29 /* move nleft, nright to right if radius decreases */
30 d1 = x[i] - x[nleft];
31 d2 = x[nright + 1] - x[i];
32 /* if d1 <= d2 with x[nright+1] == x[nright], lowest fixes */
33 if (d1 <= d2) break;
34 /* radius will not decrease by move right */
35 nleft++;
36 nright++;
38 lowest(x, y, n, x[i], &ys[i], nleft, nright, res, (iter > 1), rw, &ok);
39 /* fitted value at x[i] */
40 if (! ok) ys[i] = y[i];
41 /* all weights zero - copy over value (all rw==0) */
42 if (last < i - 1) { /* skipped points -- interpolate */
43 denom = x[i] - x[last]; /* non-zero - proof? */
44 for(j = last + 1; j < i; j = j + 1){
45 alpha = (x[j] - x[last]) / denom;
46 ys[j] = alpha * ys[i] + (1.0 - alpha) * ys[last];
49 last = i; /* last point actually estimated */
50 cut = x[last] + delta; /* x coord of close points */
51 for(i=last + 1; i < n; i++) { /* find close points */
52 if (x[i] > cut) break; /* i one beyond last pt within cut */
53 if(x[i] == x[last]) { /* exact match in x */
54 ys[i] = ys[last];
55 last = i;
58 i = max(last + 1,i - 1);
59 /* back 1 point so interpolation within delta, but always go forward */
60 } while(last < n - 1);
61 for (i = 0; i < n; i++) /* residuals */
62 res[i] = y[i] - ys[i];
63 if (iter > nsteps) break; /* compute robustness weights except last time */
64 for (i = 0; i < n; i++)
65 rw[i] = fabs(res[i]);
66 sort(rw,n);
67 m1 = 1 + n / 2; m2 = n - m1 + 1;
68 cmad = 3.0 * (rw[m1] + rw[m2]); /* 6 median abs resid */
69 c9 = .999 * cmad; c1 = .001 * cmad;
70 for (i = 0; i < n; i++) {
71 r = fabs(res[i]);
72 if(r <= c1) rw[i] = 1.0; /* near 0, avoid underflow */
73 else if(r > c9) rw[i] = 0.0; /* near 1, avoid underflow */
74 else rw[i] = pow2(1.0 - pow2(r / cmad));
77 return(0);
81 static lowest(x, y, n, xs, ys, nleft, nright, w, userw, rw, ok)
82 double *x, *y, *w, *rw, xs, *ys;
83 int n, nleft, nright, userw, *ok;
85 double range, h, h1, h9, a, b, c, r;
86 int j, nrt;
88 range = x[n - 1] - x[0];
89 h = fmax(xs - x[nleft], x[nright] - xs);
90 h9 = .999 * h;
91 h1 = .001 * h;
93 /* compute weights (pick up all ties on right) */
94 a = 0.0; /* sum of weights */
95 for(j = nleft; j < n; j++) {
96 w[j]=0.0;
97 r = fabs(x[j] - xs);
98 if (r <= h9) { /* small enough for non-zero weight */
99 if (r > h1) w[j] = pow3(1.0-pow3(r/h));
100 else w[j] = 1.0;
101 if (userw) w[j] = rw[j] * w[j];
102 a += w[j];
104 else if (x[j] > xs) break; /* get out at first zero wt on right */
106 nrt = j - 1; /* rightmost pt (may be greater than nright because of ties) */
107 if (a <= 0.0) *ok = FALSE;
108 else { /* weighted least squares */
109 *ok = TRUE;
111 /* make sum of w[j] == 1 */
112 for (j = nleft; j <= nrt; j++) w[j] = w[j] / a;
114 if (h > 0.0) { /* use linear fit */
116 /* find weighted center of x values */
117 for (j = nleft, a = 0.0; j <= nrt; j++) a += w[j] * x[j];
119 b = xs - a;
120 for (j = nleft, c = 0.0; j <= nrt; j++)
121 c += w[j] * (x[j] - a) * (x[j] - a);
123 if(sqrt(c) > .001 * range) {
124 /* points are spread out enough to compute slope */
125 b = b/c;
126 for (j = nleft; j <= nrt; j++)
127 w[j] = w[j] * (1.0 + b*(x[j] - a));
130 for (j = nleft, *ys = 0.0; j <= nrt; j++) *ys += w[j] * y[j];
134 static compar(a, b)
135 double *a, *b;
137 if (*a < *b) return(-1);
138 else if (*a > *b) return(1);
139 else return(0);
142 static sort(x, n)
143 double *x;
144 int n;
146 qsort(x, n, sizeof(double), compar);