Derive cubic spline kernel

poc/del.m

-n = 20;
+global L
+global n x y
+global rc
 
+rc = 0.2;
+L = 1;
+n = 10;
 x = rand(n, 1);
 y = rand(n, 1);
 
-T = delaunay(x, y);
-triplot(T, x, y);
+function w = wk(r)
+  global rc
+  A = 1;
+  r /= rc; r *= 2;
+  if r < 1; w = 1 - 3/2*r^2
+endfunction
+
+function r = wrap(r, r0)
+  global L
+  dr = r - r0;
+  if     abs(dr - L) < abs(dr); r -= L;
+  elseif abs(dr + L) < abs(dr); r += L; endif
+endfunction
+
+function r = dist(i, j)
+  global x y
+  xi = x(i); yi = y(i);
+  xj = x(j); yj = y(j);
+  xj = wrap(xj, xi); yj = wrap(yj, yi);
+  dx = xi - xj; dy = yi - yj;
+  r = sqrt(dx^2 + dy^2);
+endfunction
+
+d = zeros(n, 1) # density
+for i = 1:n; for j = 1:n
+  if i == j; continue; endif
+  r = dist(i, j);
+  if r > rc; continue; endif
+  w = wk(i)
+endfor; endfor
 

poc/ker/cubic.mac

+/* cubic spline kernel */
+
+v0: [ 1, b0, c0, d0];
+v1: [a1, b1, c1, d1]
+mo: [1, x, x^2, x^3];
+
+f0: v0 . mo;
+f1: v1 . mo;
+
+df0: diff(f0, x);
+df1: diff(f1, x);
+
+ddf0: diff(f0, x, 2);
+ddf1: diff(f1, x, 2);
+
+e(v, f):=subst('x = v, f);
+
+eq: [
+                    e(0,       df0),
+     e(2,      f1), e(2,       df1), e(2,        ddf1),
+     e(1, f0 - f1), e(1, df0 - df1), e(1, ddf0 - ddf1)];
+
+v : append(v0, v1);
+v : listofvars(v);
+so: linsolve(eq, v);
+
+f00: subst(so, f0);
+f10: subst(so, f1);
+
+f00: horner(f00);
+f10: factor(f10);
+
+draw2d(explicit(f00, x, 0, 1), explicit(f10, x, 1, 2)) $