PRO spline_smooth,X,Y,Yerr,distance,coefficients,smoothness,xplot,yplot, $
    TEXTOUT=textout,XTITLE=xtitle,YTITLE=ytitle, INTERP=interp,SILENT=silent, $
    PLOT=plot,ERRBAR=errbar
;+
; NAME:
;       SPLINE_SMOOTH
;
; PURPOSE:
;       Compute a cubic smoothing spline to (weighted) data
; EXPLANATION:
;       Construct cubic smoothing spline (or give regression solution) to given
;       data with minimum "roughness" (measured by the energy in the second 
;       derivatives) while restricting the weighted mean square distance
;       of the approximation from the data.  The results may be written to
;       the screen or a file or both and are optionally returned in the 
;       parameters.  The results may be optionally displayed graphically.
;
; CALLING SEQUENCE:
;       SPLINE_SMOOTH,X,Y,Yerr,distance, [coefficients,smoothness,xplot,yplot 
;               [ XTITLE= ,YTITLE=, INTERP=, TEXTOUT=,/SILENT,/PLOT,/ERRBAR]
;
; INPUT PARAMETERS:
;       X - N_POINT element vector containing the data abcissae
;       Y - N_POINT element vector containing the data ordinates
;       Yerr -     estimated uncertainty in ordinates ( positive scalar)
;       distance - upper bound on the weighted mean square distance
;                  of the approximation from the data (non-negative scalar)
;
; OPTIONAL INPUT PARAMETERS
;      xplot -    vector of spline evaluation abcissae
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
;       TEXTOUT - Controls print output device, defaults to !TEXTOUT
;       
;               textout=1       TERMINAL using /more option
;               textout=2       TERMINAL without /more option
;               textout=3       <program>.prt
;               textout=4       laser.tmp
;               textout=5       user must open file
;               textout = filename (default extension of .prt)
;
; OPTIONAL OUTPUT PARAMETERS:
;       coefficients - N_POINT x 4 element array containing the sequence of
;                      spline coefficients including the smoothed ordinates.
;       smoothness  - N_POINT element vector containing the energy in second 
;                     derivatives of approximated function.
;       yplot       - vector of evaluated spline ordinates.
;
; OPTIONAL OUTPUT KEYWORD PARAMETERS
;       /SILENT     - suppress all printing.
;       /PLOT       - display smooth curve, data ordinates and error bars
;       /ERRBAR     - display error bars
;       XTITLE      - optional title for X-axis
;       YTITLE      - optional title for Y-axis
;       INTERP      - optionally returned interpolated smooth spline
; NOTES:
;       This procedure constructs a smoothing spline according to the method
;       described in "Fundamentals of Image Processing" by A. Jain  [Prentice-
;       Hall : New Jersey 1989].
;       If the distance parameter is sufficiently large a linear regression
;       is performed, otherwise a cubic smoothing spline is constructed.
;
;       This procedure assumes regular sampling and independent identically 
;       distributed normal errors without missing data.  The data are sorted.
;       
;       SPLINE_SMOOTH uses the non-standard system variables !TEXTOUT and
;       !TEXTUNIT.
;       These can be added to one's session using the procedure ASTROLIB.
;
; COMMON BLOCKS:
;       None.
; EXAMPLE:
;       Obtain coefficients of a univariate smoothing spline fitted to data
;       X,Y assuming normally distributed errors Yerr and write the results to
;       a file.
;
;       IDL> SPLINE_SMOOTH, X, Y, Yerr, distance, coefficients, smoothness,
;            t='spline.dat'
;
;       Fit a smoothing spline to observational data.  Suppress all printing 
;       and save the smoothed ordinates in output variables. Display results.
;
;       IDL> SPLINE_SMOOTH, X, Y, Yerr, distance, coefficients, /SILENT, /PLOT
;       
; PROCEDURES CALLED:
;       Procedures TEXTOPEN, TEXTCLOSE, PLOT, PLOTERROR
;
; RESTRCTIONS:
;       This procedure is damn slow and should probably be rewritten using
;       the Cholesky decomposition.
; AUTHOR:
;       Immanuel Freedman (after A. Jain).      December, 1993
; REVISIONS
;       January 12, 1994    I. Freedman (HSTX)  Adjusted formats
;       March   14, 1994    I. Freedman (HSTX)  Improved convergence
;       March   15, 1994    I. Freedman (HSTX)  User-specified interpolates
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Call PLOTERROR instead of PLOTERR  W. Landsman      February 1999
;-

; Return on error
  ON_ERROR,2
; Constants
  FALSE = 0
  TRUE = 1
  lambda = 0.0    ; initial linear regression
  TOLERANCE = 1.0E-3
  MAXIT = 100      ; number of iterations
; Dispatch table
;
IF N_PARAMS() LT 4 THEN BEGIN
 print,'Syntax - SPLINE_SMOOTH,X,Y,Yerr,distance,[coefficients,smoothness,'
 print,'         xplot,yplot,TEXTOUT=,XTITLE=,YTITLE=,INTERP=,/SILENT,/PLOT,
 print,'         /ERRBAR]'
 RETURN
ENDIF

SZ = size(X) & NX = SZ[1]  ;Number of data points
if SZ[0] NE 1 THEN $
BEGIN
 HELP,X
 MESSAGE,'ERROR - Data matrix is not one-dimensional'
ENDIF
SZ = size(Y) & NY = SZ[1]
if SZ[0] NE 1 THEN $
BEGIN
 HELP,Y
 MESSAGE,'ERROR - Data matrix is not one-dimensional'
ENDIF
if NY NE NX THEN $
BEGIN
 HELP,Y
 MESSAGE,'ERROR - Data vector has incorrect number of elements'
ENDIF
IF distance LT 0.0 THEN MESSAGE,'ERROR - negative distance'
IF yerr     LT 0.0 THEN MESSAGE,'ERROR - negative error value'

  N_POINT = NX ; number of data points
  h = (X[N_POINT-1] - X[0])/(N_POINT -1) ;regular sampling without missing data
  IF distance LE TOLERANCE*yerr THEN distance = TOLERANCE*yerr ; interpolation
; Sort data on X
  index = sort(X) & X = X[index] & Y = Y[index]


straight: ; Linear regression [g(x)=a+bx] is solution for large distances 
; Sample averages

  average_x = TOTAL(X)/N_POINT & average_y = TOTAL(Y)/N_POINT
  average_xx = TOTAL(X*X)/N_POINT & average_xy = TOTAL(X*Y)/N_POINT
; Regression solution
  b = (average_xy - average_x*average_y)/(average_xx -average_x*average_x)
  a = average_y - b*average_x
  g = a+b*x & S = TOTAL((Y-g)*(Y-g))  ; sum of squares
  c = 0 & d = 0
  IF S LE distance THEN BEGIN
  MESSAGE,/INFORM,'Distance sufficient for linear regression'
  linear = TRUE
  smoothness = S
  GOTO,display  
  ENDIF
;
   
; The smoothing spline g(x) = a + b.dx + c.dx*dx + d.dx*dx*dx, 0 <= dx <X - X)
;                                                                        i+1 i
linear = FALSE
MESSAGE,/INFORM,'Cubic spline regression solution'
a=fltarr(N_POINT) & b=a & c=a & d=a


; Tridiagonal Toeplitz matrix

Q =fltarr(N_POINT-2,N_POINT-2) ; Jain pp. 296 ff.
Q[lindgen(N_POINT-2)*(N_POINT-1)]   = 4.0
Q[lindgen(N_POINT-3)*(N_POINT-1)+1] = 1.0
Q[lindgen(N_POINT-3)*(N_POINT-1)+N_POINT-2] = 1.0

; Lower triangular Toeplitz matrix

L = fltarr(N_POINT,N_POINT-2)
L[lindgen(N_POINT-2)*(N_POINT+1)]   = 1.0
L[lindgen(N_POINT-2)*(N_POINT+1)+1] = -2.0
L[lindgen(N_POINT-2)*(N_POINT+1)+2] = 1.0

; Define auxiliary matrices
P = (yerr*yerr) *transpose(L)#L
v = transpose(L)#y
; Iterative Newton solution of "smoothness(lambda) = S"
iteration = 0
REPEAT BEGIN
 
ainverse = P + lambda*Q 
SVDC,ainverse,svdw,svdu,svdv          ; Decompose square matrix 
small = WHERE(svdw LT TOLERANCE*MAX(svdw), count)
IF count NE 0 THEN svdw[small] = 1.0 ; Avoid division by zero below threshold
n = N_ELEMENTS(svdw) & svdwp = fltarr(n,n)
svdwp[(n+1)*lindgen(n)] = 1.0/svdw
A = svdv # svdwp # transpose(svdu)

smoothness = transpose(v)#A#P#A#v & smoothness = smoothness[0]
 derivative = 2.0 * transpose(v)#A#Q#A#P#A#v & derivative = derivative[0]
 if derivative EQ 0 THEN MESSAGE,'ERROR - convergence failed. Stop.'
 lambda = lambda + ((smoothness - distance)/derivative) ; correct sign
 OK = (ABS(smoothness - distance) LE TOLERANCE*distance)
iteration = iteration + 1
done = OK OR iteration GE MAXIT
ENDREP UNTIL done
msg = 'ERROR - Convergence failed after ' + string(MAXIT) + ' iterations.'
IF done AND NOT OK THEN BEGIN
 MESSAGE,/INFORM,strcompress(msg) 
 PRINT,'If you want to try again enter a non-negative smoothness...'
 READ,distance
 IF distance GE 0 THEN GOTO, straight
ENDIF
msg = 'Converged after ' + string(iteration) + ' iterations'
MESSAGE,/INFORM,strcompress(msg)


display:         ; evaluate and display results
IF distance LE TOLERANCE*yerr THEN distance = 0.0 ; interpolating splines
IF linear THEN BEGIN
 coefficients = fltarr(2)
 coefficients[0] = a
 coefficients[1] = b
ENDIF ELSE BEGIN

coefficients = fltarr(N_POINT, 4) 
c = lambda*A#v 
a = Y - (yerr*yerr/lambda)*L#c
c = [0.0,c,0.0]
d = (shift(c,-1) - c)/3.0*h & d[0] = 0.0
b = (shift(a,-1) -a)/h -h*c -h*h*d & b[N_POINT-1] = 0

coefficients[*, 0] = a   ; knots at each data point
coefficients[*, 1] = b
coefficients[*, 2] = c
coefficients[*, 3] = d
ENDELSE

; Open output file
IF NOT KEYWORD_SET(TEXTOUT) THEN TEXTOUT = textout
textopen,'SPLINE_SMOOTH',TEXTOUT = textout
printf,!TEXTUNIT,'SPLINE_SMOOTH: '+ SYSTIME()
; print results 
IF NOT KEYWORD_SET(silent) THEN BEGIN
; print formatted coefficients
printf,!TEXTUNIT,' '
 IF linear THEN BEGIN
  sign = " " &  IF b > 0 THEN sign = ' + '
  msg = 'Linear regression: Y = ' + string(FORMAT = '(G10.4)',a) + ' 
  msg = msg + sign + string(FORMAT='(G10.4)',b) + ' X'

 printf,!TEXTUNIT,strcompress(msg)
 ENDIF ELSE BEGIN
 msg = 'Cubic spline regression at knot points: Y = a +bX + cX^2 +d*X^3 '
 printf,!TEXTUNIT,strcompress(msg)
 msg = string(FORMAT = '(6(6X,A1,6X))','X','Y','a','b','c','d')
 printf,!TEXTUNIT,msg
  FOR i=0, N_POINT-1 DO BEGIN
    printf,!TEXTUNIT,X[i],Y[i],a[i],b[i],c[i],d[i]
  ENDFOR
  conf=[N_POINT - sqrt(2.0*N_POINT), N_POINT + sqrt(2.0*N_POINT)]
  msg = 'Confidence interval for smoothness: [' + string(conf[0]) 
  msg = msg + ', ' + string(conf[1]) + ']'
  printf,!TEXTUNIT,strcompress(msg)
 ENDELSE
 msg = 'Distance = '+string(distance) + ', Smoothness = ' + string(smoothness)
 printf,!TEXTUNIT,strcompress(msg)
ENDIF
; Close output file
textclose, TEXTOUT = textout

IF NOT KEYWORD_SET(plot) THEN BEGIN
 IF KEYWORD_SET(interp) THEN MESSAGE,/INF,'INTERP keyword ignored' 
ENDIF ELSE BEGIN

; plot results (piecewise cubic polynomial => GE 4*N_POINT+1)
 SZ = SIZE(xplot)
 sigma = replicate(yerr,N_POINT)
 IF SZ[0] EQ 0 THEN BEGIN
 MESSAGE,/INFORM,'User did not supply evaluation points'
 xplot = (0.25*h)*findgen(4*N_POINT) + x[0] ; regular sampling
ENDIF ELSE MESSAGE,/INFORM,'User supplied evaluation points'

xplot  = xplot[sort(xplot)] ; sort into increasing order
outofrange = WHERE(xplot LT min(X) OR xplot GT max(X),count)
IF count GT 0 THEN BEGIN
 msg = 'WARNING - user supplied evaluation points out of range: '
 FOR loop = 0,N_elements(outofrange) -1 DO BEGIN
  msg = msg + string(xplot[outofrange[loop]])
 ENDFOR
 MESSAGE,/INFORM,strcompress(msg)
ENDIF
inrange = WHERE((xplot GE min(X)) AND (xplot LE max(X)),count)
IF count GT 0 THEN  xplot = xplot[inrange] ; all values in range
yplot  = fltarr(N_elements(xplot))
Xindex = long((xplot-min(X))/h)  ; truncation

IF linear THEN BEGIN
 yplot = b*xplot + a
ENDIF ELSE BEGIN
 tmp = replicate(1.0,4) 
 dx  = xplot - X[xindex]                      ; fractional part
 a   = reform(tmp#a,4*N_POINT,/OVERWRITE)
 b   = reform(tmp#b,4*N_POINT,/OVERWRITE)
 c   = reform(tmp#c,4*N_POINT,/OVERWRITE)
 d   = reform(tmp#d,4*N_POINT,/OVERWRITE)
 a = a[xindex] & b = b[xindex] & c = c[xindex] & d = d[xindex]
 yplot = ((d*dx + c)*dx + b)*dx + a           ; preserve accuracy
ENDELSE
interp=fltarr(N_elements(xplot),2)
interp[*,0] = xplot
interp[*,1] = yplot
; call to PLOTERR
; Label axes as plot of Y versus X, THICK=2

IF NOT KEYWORD_SET(xtitle) THEN XTITLE='X variate'
IF NOT KEYWORD_SET(ytitle) THEN YTITLE = 'Y variate'
if linear THEN BEGIN
 rmsg = 'Linear regression' 
ENDIF ELSE BEGIN
 rmsg = 'Cubic spline regression'
ENDELSE
msg = rmsg +', distance = '+string(FORMAT = '(G10.4)',distance)+', smoothness = '+string(FORMAT = '(G10.4)',smoothness)
msg = strcompress(msg)
TITLE = msg
symbol = 7                        ; x symbol
IF N_POINT GE 20 THEN symbol = 3  ; point
IF KEYWORD_SET(errbar) THEN BEGIN
 PLOTERROR, X, Y, sigma, PSYM = symbol,xtit=xtitle,ytit=ytitle,xsty=2,ysty=2,$
        xthick=2,ythick=2,thick=2,tit=title 
ENDIF ELSE BEGIN
 PLOT, X, Y, PSYM = symbol,xtit=xitle,ytit=ytitle,xsty=2,ysty=2, $
        xthick=2,ythick=2,thick=2,tit=title
ENDELSE
OPLOT,   xplot,yplot              ; smooth curve
ENDELSE
;
RETURN
END