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