lib/kernel.c

   1 #include "xmath.h"
   2
   3 #ifndef ROOT2PI
   4 #define ROOT2PI 2.50662827463100050241
   5 #endif  ROOT2PI
   6
   7 #ifndef nil
   8 #define nil 0L
   9 #endif  nil
  10
  11 static double kernel(x, y, w, type)
  12      double x, y, w;
  13      int type;
  14 {
  15   double z, k;
  16
  17   if (w > 0.0) {
  18     z = (x - y) / w;
  19     switch (type) {
  20     case 'B':
  21       w = 2.0 * w;
  22       z = 0.5 * z;
  23       if (-0.5 < z && z < 0.5) {
  24         z = (1.0 - 4 * z * z);
  25         k = 15.0 * z * z / 8.0;
  26       }
  27       else k = 0.0;
  28       break;
  29     case 'G':
  30       w = 0.25 * w;
  31       z = 4.0 * z;
  32       k = exp(- 0.5 * z * z) / ROOT2PI;
  33       break;
  34     case 'U':
  35       w = 1.5 * w;
  36       z = .75 * z;
  37       k = (fabs(z) < 0.5) ? 1.0 : 0.0;
  38       break;
  39     case 'T':
  40       if (-1.0 <= z && z < 0.0) k = 1.0 + z;
  41       else if (0.0 <= z && z < 1.0) k = 1.0 - z;
  42       else k = 0.0;
  43       break;
  44     default:  k = 0.0; break;
  45     }
  46     k = k / w;
  47   }
  48   else k = 0.0;
  49
  50   return(k);
  51 }
  52
  53 kernel_smooth(x, y, n, width, wts, wds, xs, ys, ns, ktype)
  54      double *x, *y, width, *wts, *wds, *xs, *ys;
  55      int n, ns, ktype;
  56 {
  57   int i, j;
  58   double wsum, ysum, lwidth, lwt, xmin, xmax;
  59
  60   if (n < 1) return(1);
  61   if (width <= 0.0) {
  62     if (n < 2) return(1);
  63     for (xmin = xmax = x[0], i = 1; i < n; i++) {
  64       if (xmin > x[i]) xmin = x[i];
  65       if (xmax < x[i]) xmax = x[i];
  66     }
  67     width = (xmax - xmin) / (1 + log((double) n));
  68   }
  69
  70   for (i = 0; i < ns; i++) {
  71     for (j = 0, wsum = 0.0, ysum = 0.0; j < n; j++) {
  72       lwidth = (wds != nil) ? width * wds[j] : width;
  73       lwt = kernel(xs[i], x[j], lwidth, ktype);
  74       if (wts != nil) lwt *= wts[j];
  75       wsum += lwt;
  76       if (y != nil) ysum += lwt * y[j];
  77     }
  78     if (y != nil) ys[i] = (wsum > 0.0) ? ysum / wsum : 0.0;
  79     else ys[i] = wsum / n;
  80   }
  81   return(0);
  82 }
  83