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
  12 kernel(double x, double y, double w, int type)
  13 {
  14   double z, k;
  15
  16   if (w > 0.0) {
  17     z = (x - y) / w;
  18     switch (type) {
  19     case 'B':
  20       w = 2.0 * w;
  21       z = 0.5 * z;
  22       if (-0.5 < z && z < 0.5) {
  23         z = (1.0 - 4 * z * z);
  24         k = 15.0 * z * z / 8.0;
  25       }
  26       else k = 0.0;
  27       break;
  28     case 'G':
  29       w = 0.25 * w;
  30       z = 4.0 * z;
  31       k = exp(- 0.5 * z * z) / ROOT2PI;
  32       break;
  33     case 'U':
  34       w = 1.5 * w;
  35       z = .75 * z;
  36       k = (fabs(z) < 0.5) ? 1.0 : 0.0;
  37       break;
  38     case 'T':
  39       if (-1.0 <= z && z < 0.0) k = 1.0 + z;
  40       else if (0.0 <= z && z < 1.0) k = 1.0 - z;
  41       else k = 0.0;
  42       break;
  43     default:  k = 0.0; break;
  44     }
  45     k = k / w;
  46   }
  47   else k = 0.0;
  48
  49   return(k);
  50 }
  51
  52 int
  53 kernel_smooth(double *x, double *y, int n,
  54               double width, double *wts, double *wds,
  55               double *xs, double *ys, int ns, int 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