#|

     class20.lsp   

   - Interpolating polynomials

   - B-Splines

   - Simple smoothing spline algorithm

|#


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Newton formulation of interpolating polynomials      ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


  (defun INTERPOLATING-POLY-COEFS (x y)
    ; Returns power coefs of interpolating poly.
    (let* ((k     (length x)) 
           (a     (make-array (list k k) :initial-element 0.0))
           (coefs (repeat 0.0 k))
           )
      (dotimes (j k)               ; fill first row with y's
               (setf (aref a 0 j) (elt y j)))
      (dotimes (ii (1- k))
               (let ((i (1+ ii))  )    ; ii = i - 1
                 (dotimes (j (- k i))
                          (setf (aref a i j)
                                (/ (- (aref a ii (1+ j))
                                      (aref a ii  j))
                                   (- (elt x (+ j i))
                                      (elt x j))
                                   )))
                 ))
      (print-matrix a t :float-digits 1)
      (dotimes (i k)
               (setf (elt coefs i)
                     (aref a i 0)))
      coefs))
      
  (defun POLY-VALUE (x tau a)
     ; Evaluate poly \sum (x - t_1)...(x - t_i) a_i
     (let ((sum  (first a))
           (i    0)
           (prod 1))
       (dolist (c (rest a))
               (setf prod (* prod (- x (elt tau i))))
               (setf  sum (+ sum (* prod c)))
               (incf  i))
       sum))


   (defun INTERPOLATING-POLY (tau y)
     ; Returns interpolating poly using newton dd form.
     (let ((c (interpolating-poly-coefs tau y))  )
       #'(lambda (x)
           (poly-value x tau c))))
           

#|

  (setf tau '(1 2 3 4))                 ; interpolation grid
  
  (setf lin  (interpolating-poly-coefs tau tau))
  (setf quad (interpolating-poly-coefs tau (^ tau 2)))
  (setf cub  (interpolating-poly-coefs tau (^ tau 3)))

  (poly-value 1.5 tau lin)
  (poly-value 1.5 tau quad)
  (poly-value 1.5 tau cub)

  (let* ((y  '(1 4 -20 10 6))
         (x  (iseq (length y)))  )
    (def p (plot-function (interpolating-poly x y)
                          0 (length x)))
    (send p :add-points x y)
    )

|#




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Truncated power bases for splines                    ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


  (defun TRUNCATED-POWER-BASIS (knots x k)
    ; Power series basis for grid defined by knots evaluated at x
    ; Can do much more efficiently if the knots are equal spaced
    (bind-columns
     (outer-product x (iseq 0 k) #'^)   ; global poly terms
     (outer-product x knots
                    #'(lambda (x1 x2)
                        (if (<= x1 x2) 0
                            (^ (- x1 x2) k)))
                    )))
  
  (defun QR-OLS (x y)
    (let ((qr (qr-decomp x))  )
      (solve
           (second qr)
           (matmult (transpose (first qr)) y)
       )))

  (defun SPLINE (x y &key knots (k 3))
    ; Fits a kth order spline using the power series bases
    ; Returns function that will compute polynomial at x
    (let* ((g     (if knots knots (rest (butlast x))))
           (basis (truncated-power-basis g x k))
           (c     (qr-ols basis y))  )
      (print-matrix basis)
      (format t "y = ~{ ~5,2f~}~%" y)
      (format t "~d spline coefs = ~{ ~6,2f~}~%" (length c) c)
      #'(lambda (newx)
          (let ((x-g (mapcar #'(lambda (gi)
                                (if (<= newx gi) 0
                                    (^ (- newx gi) k)))
                             g))  )
            ; (format t "x-g = ~a~%" x-g)
            (inner-product c 
                           (append (^ newx (iseq 0 k)) x-g)
                           )))))
           
    
#|

  (print-matrix
   (truncated-power-basis '(2 6) (iseq 11) 1)   )

  (print-matrix
   (truncated-power-basis (iseq 1 9) (iseq 0 10) 1)  )

 
  (def spl (spline (iseq 10) (iseq 10)
                  :knots '( 3 6)  :k 1))
  (def p (plot-function spl 0 10))


  (let ((x (iseq 10))
        (y (uniform-rand 10))  )
    (def spl (spline (iseq 10) y :knots '(1  3 5)  :k 3))
    (def p (plot-function spl (min x) (max x)))
    (send p :add-points x y)
    )

 (let ((x (iseq 10))
        (y (uniform-rand 10))  )
    (def spl (spline (iseq 10) y :knots nil :k 3))
    (def p (plot-function spl (min x) (max x)))
    (send p :add-points x y)
    )


|#



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    B-splines                                            ;;
;;                     deBoor, p131                        ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


  (defun B-SPLINE (x i k tau)
    ; B_{i,k}(x) using grid defined by tau
    (labels ((bsp (x i k)
                  (if (= 1 k)
                      (if (<= (elt tau i) x (elt tau (1+ i)))
                          1 0)
                      (+ (* (/ (- x (elt tau i))
                               (- (elt tau (+ i k -1))
                                  (elt tau i)))
                            (bsp x i (1- k)))
                         (* (/ (- (elt tau (+ i k)) x)
                               (- (elt tau (+ i k))
                                  (elt tau (+ i 1))))
                            (bsp x (1+ i) (1- k)))
                         )))  )
      (bsp x i k)
      ))
                                  
                                
       
#|

  (def grid (iseq 10))
  (def k 4)

  (def p
       (plot-function #'(lambda (x) (b-spline x 3 k grid))
                      (min grid) (max grid)))
       
  (let ((x (rseq (min grid) (max grid) 50))  )
    (send p :add-lines
          x
          (mapcar #'(lambda (x) (b-spline x 5 k grid)) x)
          :color 'red))
                    

|#


|#