; file: splinapp_set.pro = define a cubic spline approximation ; init: Apr 18 2016 Rob Rutten Oslo = start from splinapp_nv.pro ; last: Apr 18 2016 Rob Rutten Oslo ;+ pro splinapp_set,x,y,a,b,c,d,smooth=smooth ; PURPOSE: from Nikola Vitas splinesmooth: This procedure computes ; coefficients of approximating cubic splines for a given ; observational set and smoothing parameter. The method is coded ; according to Pollock D.S.G. (1999), "A Handbook of Time-Series ; Analysis, Signal Processing and Dynamics, Academic Press", San ; Diego (60$ at Amazon on Apr 18 2016, RR). RR: Pollock's ; "Smoothing with Cubic Splines" is freely downloadable at ; ResearchGate. ; ; CATEGORY: ; OAP ; ; CALLING SEQUENCE: ; see above ; ; INPUTS: ; x,y: 1D input vectors ; ; OPTIONAL KEYWORD INPUTS: ; smooth: smooth parameter, 0 = interpolation, large = linear (default 1) ; ; OUTPUT: ; a,b,c,d: 1D arrays same size as input with spline coefficients ; in the input nodes. ; ; MODIFICATION HISTORY: ; Apr 18 2016 RR: revamp splinapp_nv.pro from Nikola Vitas ;- ; defaults for keywords if (n_elements(smooth) eq 0) then smooth=1.0 NUM = SIZE(X, /N_ELEMENTS) N = NUM-1 ; define help arrays H = DBLARR(NUM) & R = DBLARR(NUM) & F = DBLARR(NUM) & P = DBLARR(NUM) Q = DBLARR(NUM) & U = DBLARR(NUM) & V = DBLARR(NUM) & W = DBLARR(NUM) ; define output arrays A = DBLARR(NUM) & B = DBLARR(NUM) & C = DBLARR(NUM) & D = DBLARR(NUM) ; compute start values H[0] = X[1] - X[0] R[0] = 3.D/H[0] ; compute other H, R, F, P & Q FOR I = 1, N - 1 DO BEGIN H[I] = X[I+1] - X[I] R[I] = 3.D/H[I] F[I] = -(R[I-1] + R[I]) P[I] = 2.D * (X[I+1] - X[I-1]) Q[I] = 3.D * (Y[I+1] - Y[I])/H[I] - 3.D * (Y[I] - Y[I-1])/H[I-1] ENDFOR ; get diagonals of matrix W + LAMBDA T' T FOR I = 1, N - 1 DO BEGIN U[I] = R[I-1]^2 + F[I]^2 + R[I]^2 U[I] = smooth * U[I] + P[I] V[I] = F[I] * R[I] + R[I] * F[I+1] V[I] = smooth * V[I] + H[I] W[I] = smooth * R[I] * R[I+1] ENDFOR ; decompose in form L' D L V[1] = V[1]/U[1] W[1] = W[1]/U[1] FOR J = 2, N-1 DO BEGIN U[J] = U[J] - U[J-2] * W[J-2]^2 - U[J-1] * V[J-1]^2 V[J] = (V[J] - U[J-1] * V[J-1] * W[J-1])/U[J] W[J] = W[J]/U[J] ENDFOR ; Gaussian elimination to solve Lx = T'y Q[0] = 0.D FOR J = 2, N-1 DO Q[J] = Q[J] - V[J-1] * Q[J-1] - W[J-2] * Q[J-2] FOR J = 1, N-1 DO Q[J] = Q[J]/U[J] ; Gaussian elimination to solve L'c = D^{-1}x Q[N-2] = Q[N-2] - V[N-2]*Q[N-1] FOR J = N-3, 1, -1 DO Q[J] = Q[J] - V[J] * Q[J+1] - W[J] * Q[J+2] ; coefficients for the first segment ;RR not so good D[0] = Y[0] - smooth * R[0] * Q[1] D[1] = Y[1] - smooth * (F[1] * Q[1] + R[1] * Q[2]) A[0] = Q[1]/(3.D * H[0]) B[0] = 0.D C[0] = (D[1] - D[0])/H[0] - Q[1] * H[0]/3.D ; other coefficient FOR J = 1, N-1 DO BEGIN A[J] = (Q[J+1]-Q[J])/(3.D * H[J]) B[J] = Q[J] C[J] = (Q[J] + Q[J-1]) * H[J-1] + C[J-1] D[J] = R[J-1] * Q[J-1] + F[J] * Q[J] + R[J] * Q[J+1] D[J] = Y[J] - smooth * D[J] ENDFOR D[N] = Y[N] - smooth * R[N-1] * Q[N-1] END ; ================= test per IDLWAVE Hyper-C ========================= ; input: Shanghai rankings for Utrecht "University" against time year=[2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015] rank=[ 40, 39, 41, 40, 42, 47, 52, 50, 48, 53, 52, 57, 56] year=year+8./12 ; announcements come on August 15 nyr=n_elements(year) splinapp_set,year,rank,a,b,c,d ;; ; get output at input abscissae ;; forward_function splinapp_get ;; splinfit=rank ;; for ix=0,nyr-1 do splinfit[ix]=splinapp_get(year[ix],year,a,b,c,d) ;; ; make plot ;; xrange=[fix(year[0]),fix(year[nyr-1])+2] ;; yrange=[60,33] ;; plot,year,rank,psym=6,$ ;; position=[0.2,0.2,0.8,0.7],/normal,$ ; set margins around plot ;; xticklen=0.03,yticklen=0.03/1.6,$ ; same-length ticks ;; xtitle=xtitle,ytitle=ytitle,$ ;; xrange=xrange,xstyle=1,yrange=yrange,ystyle=1 ;; ; overplot cubic spline approximation ;; oplot,year,splinfit ; same for many samples tryx=indgen(100)*16./100+2002 splinfit=fltarr(100) forward_function splinapp_get for ix=0,99 do splinfit[ix]=splinapp_get(tryx[ix],year,a,b,c,d) ; make plot xrange=[2002,2018] yrange=[60,33] plot,year,rank,psym=6,$ position=[0.2,0.2,0.8,0.7],/normal,$ ; set margins around plot xticklen=0.03,yticklen=0.03/1.6,$ ; same-length ticks xtitle=xtitle,ytitle=ytitle,$ xrange=xrange,xstyle=1,yrange=yrange,ystyle=1 ; overplot cubic spline approximation oplot,tryx,splinfit end