### Thread: 3 splines in j

1. #### 3 splines in j

The j programming language may be installed from J Home .
Algorithm for leastSquare is converted to a locale (class)

Of the graphs (see end of post), note the extrapolation of the cosine curves for the natural and notaknot curves which demonstrate the endpoint conditions. Extrapolation is not recommended.

spline.ijs
Code:
```coclass'natural'

doc=: 0 :0

spline generates spline coefficients. x are independent, y are dependent
Since this is a class generate the spline object and evaluate it with

instance =: (x ,: y) conew 'spline'
eval__instance independent_value

This makes a cubic spline matching value and slope at the knots,
but for the end points which have zero concavity.

https://en.wikipedia.org/wiki/Spline_interpolation
)

create=: 3 :0
'KNOTS Y'=: y
assert # KNOTS
assert (-: ~.) KNOTS  NB. knots are unique
assert (-: /:~) KNOTS NB. knots are ordered for interval index
dx =. 2 -~/\ KNOTS
rdx =: % dx
dy =. 2 -~/\ Y
A =. (1 2 1 *"1 (}:  ,. (}: + }.) ,. }.)) rdx
A =. (-@:i.@:# |."0 1 ({.~ (0 2 + #))) A
A =. (2 1 * {. rdx) , A , (- 2 + # A) {. 1 2 * {: rdx
c =. dy * *: rdx
b =: (3 * 0&, + ,&0) c
k =: b %. A
a =: ((}: k) * dx) - dy
b =: dy - (}. k) * dx
)

eval=: 3 :0&>
i =: (y = {. KNOTS) + KNOTS I. y
if. i e. 0 , # KNOTS do.
if. i = 0 do. f =. {. else. f =. {: end.
slope =. f k
why =. f Y
ex =. f KNOTS
why + slope * y - ex
return.
end.
i =: <: i
t =: (y - i { KNOTS) * i { rdx
w =: -. t
((Y {~ 0 1 + i) +/ .* w , t) + t * w * (((i{a) * w) + (i{b) * t)
)

destroy=: codestroy@:(''"_)

EXAMPLE=: noun define
require'plot'
cos =: 2&o.
O =: conew&'natural'(,: cos) 1p1 * i: 1j20
'output cairo 600 300'plot (j. eval__O)  4 * i: 1j100
)

coclass'notaknot'

doc=: 0 :0

spline generates spline coefficients. x are independent, y are dependent
Since this is a class generate the spline object and evaluate it with

instance =: (x ,: y) conew 'spline'
eval__instance independent_value

This makes a cubic spline matching value and slope at the knots,
but for the end points which have zero third derivative.

https://en.wikipedia.org/wiki/Spline_interpolation
)

create=: 3 :0
'KNOTS Y'=: y
assert # KNOTS
assert (-: ~.) KNOTS  NB. knots are unique
assert (-: /:~) KNOTS NB. knots are ordered for interval index
dx =. 2 -~/\ KNOTS
rdx =: % dx
dy =. 2 -~/\ Y
A =. (1 2 1 *"1 (}:  ,. (}: + }.) ,. }.)) rdx
A =. (-@:i.@:# |."0 1 ({.~ (0 2 + #))) A
A =. ((*:0{rdx),((+*-)/2{.rdx),-*:1{rdx),A,(- 2 + # A) {. (*:_2{rdx),((+*-)/_2{.rdx),-*:{:rdx
c =. dy * *: rdx
b =: (3 * 0&, + ,&0) c
head =. +: ((0{dy)*3^~0{rdx) - ((1{dy)*3^~1{rdx)
tail =. +: ((_2{dy)*3^~_2{rdx) - ((_1{dy)*3^~_1{rdx)
b =:  _1 tail"_`[`]} b
k =: b %. A
a =: ((}: k) * dx) - dy
b =: dy - (}. k) * dx
)

eval=: 3 :0&>
i =: (y = {. KNOTS) + KNOTS I. y
if. i e. 0 , # KNOTS do.
if. i = 0 do. f =. {. else. f =. {: end.
slope =. f k
why =. f Y
ex =. f KNOTS
why + slope * y - ex
return.
end.
i =: <: i
t =: (y - i { KNOTS) * i { rdx
w =: -. t
((Y {~ 0 1 + i) +/ .* w , t) + t * w * (((i{a) * w) + (i{b) * t)
)

destroy=: codestroy@:(''"_)

EXAMPLE=: noun define
require'plot'
cos =: 2&o.
O =: conew&'notaknot'(,: cos) 1p1 * i: 1j20
'output cairo 600 300'plot (j. eval__O)  4 * i: 1j100
)

coclass'leastSquare'

doc=: 0 :0
https://code.jsoftware.com/mediawiki/images/c/cb/Ls_spline.ijs

fit spline to noisy data

locale leastSquare   find cubic spline fit by least squares.
From http://www.astro.umd.edu/~jph/ls_spline.ijs by
J. Patrick Harrington (jph@astro.umd.edu).
Ref.: S K Lucas, Gazette Aust. Math. Soc., 30, 207 (2003).
t= boundaries of cubic segments, (xi,yi) = data points.
Use: OBJECT=: conew_leastSquare_ t;xi;yi
eval__OBJECT locations
)

x =. > {. y
y =. > }. y
'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
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
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
zp=: (Np1+i. Np1){ZZ
>t; z; zp
)

NB. Thanks to Arie Groeneveld for improvements 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=: [: <: (1"_ >. [: +/"1 ] >/ [) <. [: <: [: # [

NB. Least squares spline evaluation using object
NB. of leastSquare  eval__OBJECT y --> f(y)
NB. The argument y may be an array of x-values.

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
)

destroy=: codestroy@:(''"_)

NB. plot dots & line: (dots x;y) dot_and_line (line x;y)
require'plot'
pd 'output cairo 600 300' [ pd'reset'     NB. Size 300x600 pixels
pd'type dot;pensize 4;color red'          NB. Distinguish dots as red
pd x
pd'type line;pensize 3; color blue'
pd y
pd'show'
)

NB. Test using classic "titanium heat data":
EXAMPLE=: noun define
int01=: i.@>:%]  NB. int01 n --> divides [0 ... 1] into n intervals.

ti=:    0.644 0.622 0.638 0.649 0.652 0.639 0.646 0.657 0.652 0.655
ti=: ti,0.644 0.663 0.663 0.668 0.676 0.676 0.686 0.679 0.678 0.683
ti=: ti,0.694 0.699 0.710 0.730 0.763 0.812 0.907 1.044 1.336 1.881
ti=: ti,2.169 2.075 1.598 1.211 0.916 0.746 0.672 0.627 0.615 0.607
ti=: ti,0.606 0.609 0.603 0.601 0.603 0.601 0.611 0.601 0.608

xx=. 595+ 10*i.#ti
NB. 8 interval points "judiciously" chosen (for best fit!):
t=. 595 820 850 870 900 915 980 1075
O =: conew&'leastSquare' t;xx;ti
uu=. 595+(1075-595)*int01 400
vv=. eval__O uu
(xx;ti) dot_and_line__O uu;vv  NB. this is the plot shown as "Fit ..."
)```
natural.png
notaknot.png
leastSquare.png
Last edited by b49P23TIvg; March 16th, 2019 at 06:17 PM.