lib/lowess.c

   1 /*
   2   Translated from RATFOR lowess code of W. S. Cleveland as obtained from NETLIB
   3 */
   4
   5 #include "xmath.h"
   6 /*
   7 #ifdef __APPLE__
   8 #include <stdlib.h>
   9 #endif
  10 */
  11 #include <stdlib.h>
  12
  13 #define FALSE 0
  14 #define TRUE 1
  15
  16 extern long max(long, long);
  17 extern long min(long, long);
  18
  19 static double pow2(double x) { return(x * x); }
  20 static double pow3(double x) { return(x * x * x); }
  21 /* static */ double fmax(double x, double y) { return (x > y ? x : y); }
  22
  23 int
  24 static compar(const void *aa, const void *bb)
  25 {
  26   const double* a = (double*)aa;
  27   const double* b = (double*)bb;
  28
  29   if (*a < *b) return(-1);
  30   else if (*a > *b) return(1);
  31   else return(0);
  32 }
  33
  34 static void
  35 lowest(double *x, double *y, size_t n, double xs, double *ys, long nleft, long nright,
  36        double *w, int userw, double *rw, int *ok)
  37 {
  38   double range, h, h1, h9, a, b, c, r;
  39   long j, nrt;
  40
  41   range = x[n - 1] - x[0];
  42   h = fmax(xs - x[nleft], x[nright] - xs);
  43   h9 = .999 * h;
  44   h1 = .001 * h;
  45
  46   /* compute weights (pick up all ties on right) */
  47   a = 0.0;        /* sum of weights */
  48   for(j = nleft; j < n; j++) {
  49     w[j]=0.0;
  50     r = fabs(x[j] - xs);
  51     if (r <= h9) {    /* small enough for non-zero weight */
  52       if (r > h1) w[j] = pow3(1.0-pow3(r/h));
  53       else w[j] = 1.0;
  54       if (userw) w[j] = rw[j] * w[j];
  55       a += w[j];
  56     }
  57     else if (x[j] > xs) break;  /* get out at first zero wt on right */
  58   }
  59   nrt = j - 1;  /* rightmost pt (may be greater than nright because of ties) */
  60   if (a <= 0.0) *ok = FALSE;
  61   else { /* weighted least squares */
  62     *ok = TRUE;
  63
  64     /* make sum of w[j] == 1 */
  65     for (j = nleft; j <= nrt; j++) w[j] = w[j] / a;
  66
  67     if (h > 0.0) {     /* use linear fit */
  68
  69       /* find weighted center of x values */
  70       for (j = nleft, a = 0.0; j <= nrt; j++) a += w[j] * x[j];
  71
  72       b = xs - a;
  73       for (j = nleft, c = 0.0; j <= nrt; j++)
  74         c += w[j] * (x[j] - a) * (x[j] - a);
  75
  76       if(sqrt(c) > .001 * range) {
  77         /* points are spread out enough to compute slope */
  78         b = b/c;
  79         for (j = nleft; j <= nrt; j++)
  80           w[j] = w[j] * (1.0 + b*(x[j] - a));
  81       }
  82     }
  83     for (j = nleft, *ys = 0.0; j <= nrt; j++) *ys += w[j] * y[j];
  84   }
  85 }
  86
  87 static void
  88 sort(double *x, size_t n)
  89 {
  90   extern void qsort();
  91
  92   qsort(x, n, sizeof(double), compar);
  93 }
  94
  95
  96
  97 int
  98 lowess(double *x, double *y, size_t n,
  99        double f, size_t nsteps,
 100        double delta, double *ys, double *rw, double *res)
 101 {
 102   int iter, ok;
 103   long i, j, last, m1, m2, nleft, nright, ns;
 104   double d1, d2, denom, alpha, cut, cmad, c9, c1, r;
 105
 106   if (n < 2) { ys[0] = y[0]; return(1); }
 107   ns = max(min((long) (f * n), n), 2);  /* at least two, at most n points */
 108   for(iter = 1; iter <= nsteps + 1; iter++){      /* robustness iterations */
 109     nleft = 0; nright = ns - 1;
 110     last = -1;        /* index of prev estimated point */
 111     i = 0;   /* index of current point */
 112     do {
 113       while(nright < n - 1){
 114         /* move nleft, nright to right if radius decreases */
 115         d1 = x[i] - x[nleft];
 116         d2 = x[nright + 1] - x[i];
 117         /* if d1 <= d2 with x[nright+1] == x[nright], lowest fixes */
 118         if (d1 <= d2) break;
 119         /* radius will not decrease by move right */
 120         nleft++;
 121         nright++;
 122       }
 123       lowest(x, y,
 124              n, x[i],
 125              &ys[i],
 126              nleft, nright,
 127              res, (iter > 1), rw, &ok);
 128       /* fitted value at x[i] */
 129       if (! ok) ys[i] = y[i];
 130       /* all weights zero - copy over value (all rw==0) */
 131       if (last < i - 1) { /* skipped points -- interpolate */
 132         denom = x[i] - x[last];    /* non-zero - proof? */
 133         for(j = last + 1; j < i; j = j + 1){
 134           alpha = (x[j] - x[last]) / denom;
 135           ys[j] = alpha * ys[i] + (1.0 - alpha) * ys[last];
 136         }
 137       }
 138       last = i;        /* last point actually estimated */
 139       cut = x[last] + delta;     /* x coord of close points */
 140       for(i=last + 1; i < n; i++) {     /* find close points */
 141         if (x[i] > cut) break;     /* i one beyond last pt within cut */
 142         if(x[i] == x[last]) {      /* exact match in x */
 143           ys[i] = ys[last];
 144           last = i;
 145         }
 146       }
 147       i = max(last + 1,i - 1);
 148       /* back 1 point so interpolation within delta, but always go forward */
 149     } while(last < n - 1);
 150     for (i = 0; i < n; i++)      /* residuals */
 151       res[i] = y[i] - ys[i];
 152     if (iter > nsteps) break; /* compute robustness weights except last time */
 153     for (i = 0; i < n; i++)
 154       rw[i] = fabs(res[i]);
 155     sort(rw,n);
 156     m1 = 1 + n / 2; m2 = n - m1 + 1;
 157     cmad = 3.0 * (rw[m1] + rw[m2]);      /* 6 median abs resid */
 158     c9 = .999 * cmad; c1 = .001 * cmad;
 159     for (i = 0; i < n; i++) {
 160       r = fabs(res[i]);
 161       if(r <= c1) rw[i] = 1.0;      /* near 0, avoid underflow */
 162       else if(r > c9) rw[i] = 0.0;  /* near 1, avoid underflow */
 163       else rw[i] = pow2(1.0 - pow2(r / cmad));
 164     }
 165   }
 166   return(0);
 167 }
 168
 169
 170