NB. Routine to find cubic spline fit to data by least squares.
NB. Ref.: S K Lucas, Gazette Aust. Math. Soc., 30, 207 (2003).   
NB. t= boundaries of cubic segments, (xi,yi) = data points.
NB. Usage: 't z zp'=. OUT=. t ls_spline xi;yi
NB. z and zp are values of spline fit f(x) and df/dx at t's.
NB. (z & zp are sufficient to evaluate the fit f(x) at any x.)

ls_spline=: dyad define 
'xi yi'=. /:~&.|: >y   NB. sort xi into an ascending sequence 
idx=. x IdX xi         NB. idx is t interval in which each xi resides
NB. if first or last interval(s) are empty of data, delete them
t=. ({. idx)}.((2-#x)+{:idx)}. x 
Np1=. 2+ Nm1=. <: N=. #h=. (}.-}:)t
bx=. <"1 |: >(i. n=. #xi); t IdX xi
v=. _1+ u=. bx{ h%"1~ xi -/ }:t
uu=. *: u=. u bx}0$~n,N
vv=. *: v=. v bx}0$~n,N
aa=. +/ *:al=. (1+2*u)*vv
ab=. +/ al*bt=. uu*(1-2*v)
ag=. +/ al*gm=. h*"1 u*vv
ad=. +/ al*dl=. h*"1 uu*v 
bd=. +/ bt*dl [bg=. +/ bt*gm [bb=. +/ *:bt
dd=. +/ *:dl  [gd=. +/ gm*dl [gg=. +/ *:gm
D=. +:(0,~ +/ yi*al)+(0, +/ yi*bt)
G=. +: (0,~ +/ yi*gm)+(0, +/ yi*dl)
bx11=. >: &.> bx00=. ;/ ,.~i. N        
bx01=. (0 1)&+ &.> bx00
bx10=. (1 0)&+ &.> bx00
A=. (bb bx11}E) + aa bx00}E=. 0$~,~Np1
A=. +: (ab,ab) (bx01,bx10)}A
B=. (ag bx00}E) + bd bx11}E
B=. +: (bg,ad) (bx01,bx10)}B
E=. (gg bx00}E) + dd bx11}E
E=. +: (gd,gd) (bx01,bx10)}E
bx00=. }: bx00 [ bx01=. }: bx01
bx02=. (0 2)&+ &.> bx00
hip=. }.h  [hi=. }:h 
C=. (3*hip%hi) bx00}F=. 0$~Nm1,Np1
C=. (3*(hi%hip)-hip%hi) bx01}C
C=. (_3*hi%hip) bx02}C
F=. hip bx00}F
F=. (+: hi+hip) bx01}F
F=. hi bx02}F
NB. set up matrix of linear system to be solved:
M=. (A,.(|:B),.|:C),(B,.E,.|:F),C,.F 
RHS=. (D,G,Nm1#0) NB. right hand side of system
ZZ=. RHS %. M     NB. solve for f(x_i) & [df/dx]_i
z=. (i. Np1){ZZ
>t; z; zp=. (Np1+i. Np1){ZZ
)

NB. Thanks to Arie Groeneveld for imporvements to this code.
NB. Note: equation (2) of Lucas should be u_i(x)=(x-t_i)/h_i
NB. and the 2nd B in the pseudocode matrix should not be B^T
NB. but just B as in equation (7).

NB. t IdX x0 --> index of where x0 fits in t-array:
NB.    0   1   2   3   4   5   6  <- t values
NB.    |   |   |   |   |   |   |
NB.  <-- 0 | 1 | 2 | 3 | 4 | 5 --->  IdX values
IdX=: ] I.~ [:}. [:}: [    NB. (}.}: t) I. x0

NB. Least squares spline evaluation using output "OUT"
NB. of ls_spline. Usage: OUT eval_spline y --> f(y)
NB. The argument y may be an array of x-values.

eval_spline=: dyad define
't z zp'=. x
NB. Array t holds the boundaries of the spline pieces, z is
NB. the array of the values of the cubic splines at the
NB. t points, and zp (z') the derivatives at those points.
N=. # h=. (}.-}:) t
v=. _1+ u=. h%~ |: y -/ }:t
a=. (1+2*u)*v*v* }:z 
b=. u*u*(1-2*v)* }.z 
c=. h*u*v*v* }:zp 
d=. h*u*u*v* }.zp 
box=. <"1 |: >(i. #y); t IdX y
box { ((#y),N)$,|: a+b+c+d
)

NB. plot dots & line: (dots x;y) dot_and_line (line x;y)
require'plot'
dot_and_line=: dyad define
pd'reset'
pd'type dot;pensize 4'
pd x
pd'type line;pensize 3'
pd y
pd'show'
)

NB. Example:
NB. t=. 0, 0.2 0.5 0.7 1
NB. xx=. rand 35
NB. yy=. (cos 1p1*xx)+_0.05+0.1*rand 35
NB. 'type symbol;pensize 3' plot xx;yy
NB. load'ls_spline.ijs'
NB. ]OUT=. t ls_spline xx;yy
NB. uu=. int01 300
NB. vv=. OUT eval_spline uu
NB. (xx;yy) dot_and_line uu;vv

rand=: [: ? ] # 0"_
int01=: i.@>:%]  NB. int01 n --> divides [0 ... 1] into n intervals.

NB. Test using classic "titanium heat data":
NB. ti=:    0.644 0.622 0.638 0.649 0.652 0.639 0.646 0.657 0.652 0.655
NB. ti=: ti,0.644 0.663 0.663 0.668 0.676 0.676 0.686 0.679 0.678 0.683
NB. ti=: ti,0.694 0.699 0.710 0.730 0.763 0.812 0.907 1.044 1.336 1.881
NB. ti=: ti,2.169 2.075 1.598 1.211 0.916 0.746 0.672 0.627 0.615 0.607
NB. ti=: ti,0.606 0.609 0.603 0.601 0.603 0.601 0.611 0.601 0.608
NB.
NB. xx=. 595+ 10*i.#ti
NB. 8 interval points "judiciously" chosen (for best fit!):
NB. t=. 595 820 850 870 900 915 980 1075
NB. OUT=. t ls_spline xx;ti
NB. uu=. 595+(1075-595)*int01 400
NB. vv=. OUT eval_spline uu
NB. (xx;ti) dot_and_line uu;vv  NB. this is the plot shown as "Fit ..."

========================================================================

NB. Least squares fit of points to "not quite a cubic spline".
NB. see Ferguson and Staley, Comm. ACM 16, 380 (1973).
NB. This routine is simpler than ls_spline above, since it doesn't enforce
NB. continuity of the second derivatives. The results usually look similar 
NB. but are not so smooth. Given for comparison, but I don't recommend it. 
NB. t= frets, xi and yi data points. 
NB. Usage: 't z zp'=. OUT2=. t ls_spline2 xi;yi
NB. z and zp are values of spline fit f(x) and df/dx at t's.

ls_spline2=: dyad define 
'xi yi'=. /:~&.|: >y   NB. sort xi into an ascending sequence 
idx=. x IdX xi
NB. if first or last interval(s) are empty of data, delete them
t=. ({. idx)}.((2-#x)+{:idx)}. x 
Np1=. >: N=. #h=. (}.-}:)t
bx=. <"1 |: >(i. n=. #xi); t IdX xi
v=. _1+ u=. bx{ h%"1~ xi -/ }:t
uu=. *: u=. u bx}0$~n,N
vv=. *: v=. v bx}0$~n,N
aa=. +/ *:al=. (1+2*u)*vv
ab=. +/ al*bt=. uu*(1-2*v)
ag=. +/ al*gm=. h*"1 u*vv
ad=. +/ al*dl=. h*"1 uu*v
bd=. +/ bt*dl [bg=. +/ bt*gm [bb=. +/ *:bt
dd=. +/ *:dl  [gd=. +/ gm*dl [gg=. +/ *:gm
D=. +:(0,~ +/ yi*al)+(0, +/ yi*bt)
G=. +: (0,~ +/ yi*gm)+(0, +/ yi*dl)
bx11=. >: &.> bx00=. ;/ ,.~i. N        
bx01=. (0 1)&+ &.> bx00
bx10=. (1 0)&+ &.> bx00
A=. (bb bx11}E) + aa bx00}E=. 0$~,~Np1
A=. +: (ab,ab) (bx01,bx10)}A
B=. (ag bx00}E) + bd bx11}E
B=. +: (bg,ad) (bx01,bx10)}B
E=. (gg bx00}E) + dd bx11}E
E=. +: (gd,gd) (bx01,bx10)}E
M=. (A,.|:B),(B,.E)
ZZ=. (D,G)%. M  NB. solve for f(x_i) & [df/dx]_i
z=. (i. Np1){ZZ
>t; z; zp=. (Np1+i. Np1){ZZ
)

NB. Example:
NB. t=. 0, 0.2 0.5 0.7 1
NB. xx=. rand 35
NB. yy=. (cos 1.5p1*xx)+_0.05+0.1*rand 35
NB. 'type symbol;pensize 3' plot xx;yy
NB. load'ls_spline2.ijs'
NB. ]OUT2=. t ls_spline2 xx;yy
NB. uu=. int01 300
NB. vv2=. OUT2 eval_spline uu
NB. (xx;yy) dot_and_line uu;vv2
NB. Compare to ls_spline from above:
NB. (xx;yy) dot_and_line uu;>vv;vv2