clean up of consistent-dataframe-p generic/methods. They were subtly wrong.
[CommonLispStat.git] / lib / kernel.c
blob5504a47e3b7dc3dcf346443327a007e7b9b881b8
1 #include "xmath.h"
3 #ifndef ROOT2PI
4 #define ROOT2PI 2.50662827463100050241
5 #endif /* ROOT2PI*/
7 #ifndef nil
8 #define nil 0L
9 #endif /*nil */
12 #ifdef __APPLE__
13 #include <stdlib.h>
14 #endif
17 #include <stdlib.h>
19 static double
20 kernel(double x, double y, double w, int type)
22 double z, k;
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;
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;
53 k = k / w;
55 else k = 0.0;
57 return(k);
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)
65 size_t i, j;
66 double wsum, ysum, lwidth, lwt, xmin, xmax;
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];
75 width = (xmax - xmin) / (1 + log((double) n));
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];
86 if (y != nil) ys[i] = (wsum > 0.0) ? ysum / wsum : 0.0;
87 else ys[i] = wsum / n;
89 return(0);