Thread: 3 splines in j

    #1
  1. Contributing User
    Devshed God 1st Plane (5500 - 5999 posts)

    Join Date
    Aug 2011
    Posts
    5,988
    Rep Power
    510

    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 =:   0 head"_`[`]} b
     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
    )
    
    create=: monad define
     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
     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
     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.
    
    eval=: monad define
     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'
    dot_and_line=: dyad define
     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.
    [code]Code tags[/code] are essential for python code and Makefiles!

IMN logo majestic logo threadwatch logo seochat tools logo