#|

     class13.lsp   

    - Move to version 3.48 has fixed slider error

    - Regression and sweep operations
    
    - Matrix factorizations


|#

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Slider animation                                     ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;--- Function to generate test data  (saves truth as side effect)

   (defun SINE-GEN (n sigma)
     (+ (setf TRUTH (sin (rseq 0 (* 2 pi) n)))
        (* sigma (normal-rand n))
        ))


;--- Pick parameters

(def n 100)
(def noiseSD .5)


;--- Plot data and true function to estimate

(def xax (iseq n))
(def  thePlot (plot-points xax (setf DATA (sine-gen n noiseSD)) ))
(send thePlot :add-lines xax truth :color 'magenta)


;--- Moving average estimator, revised.

   (defun MOVING-AVERAGE (x width)
     "Smooths assuming input x is equally spaced."
     (unless (oddp width)
             (format t "Width needs to be odd integer >= 3.~%")
             (setf width (let ((rw (round width))  )
                           (if (oddp rw) rw (1+ rw))))   )
     (let* ((half  (/ (1- width) 2))             ; invariants done once
            (n     (length x))
            (count half)
            (sum   (sum (select x (iseq half)))) ; initialize terms
            (drop  (- (1+ half)))
            (add   half)
            (smth  ())  )
       (dotimes (i n)
                (when (<= 0 drop)
                      (decf sum (select x drop))
                      (decf count))
                (when (<  add n)
                      (incf sum (select x add))
                      (incf count))
                (incf drop)
                (incf add)
                ; (format t "~d ~d  count=~d~%" drop add count)
                (push (/ sum count) smth))
       (nreverse smth)))


(compile 'moving-average)


(interval-slider-dialog            
 (list 5 n)      ; limits on interval
 :points 50      ; divisions of interval
 :action (lambda (x)
           (send thePlot :clear-lines)
           (send thePlot :add-lines xax (moving-average data x))
           ))
         

;--- Use the C function

    (defun MOV-AVG (x width)
      (unless (oddp width)
              (format t "Width needs to be odd integer >= 3.~%")
              (setf width (let ((rw (round width))  )
                            (if (oddp rw) rw (1+ rw))))   )
      (third (call-cfun "MovAvg"
                 (float x) (length x)
                 width
                 (repeat 7.7 (length x))
                 ))
      )



(interval-slider-dialog            
 (list 5 n)      ; limits on interval
 :points 50      ; divisions of interval
 :action (lambda (x)
           (send thePlot :clear-lines)
           (send thePlot :add-lines xax (mov-avg data x))
           ))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Sweep operations and regression                      ;;
;;                              Tierney, Section 5.5, 170  ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Utility for building matrix of polynomials

   (defun POLY-MATRIX (x order)
     "Builds matrix of columns x,x^2,...."
     (apply
      #'bind-columns
      (mapcar #'(lambda (p) (^ x p)) (iseq 1 order))
      ))
 
(print-matrix 
 (setf poly (poly-matrix (rseq -1 1 100) 4)))


;--- Utilities for extracting a row or column as a list

    (defun AS-LIST (x)
      (coerce (compound-data-seq x) 'list))
 
    (defun COL-AS-LIST (matrix j)
      (as-list 
       (select matrix (iseq (array-dimension matrix 0)) j)
       ))

    (defun ROW-AS-LIST (matrix i)
      (as-list 
       (select matrix i (iseq (array-dimension matrix 1)))
       ))


(plot-points (col-as-list poly 1) (col-as-list poly 3))


;--- Make a sweep matrix "sweep"

    (setf haveSwept ())    

    (defun SWEEP (j)
      (let ((new (sweep-operator sweep (list j)))   )
        (format t "Swept rows ~a~%" (second new))
        (setf haveSwept
              (when (second new)
                    (if (member j haveSwept)
                        (remove-if #'(lambda (x) (= x j)) haveSwept)
                        (push j haveSwept))  ))
        (setf sweep (first new))
        ))

    (defun COEFS ()
      (format t "Coefs = ~a ~% ~a ~%" haveSwept
               (select
                (row-as-list sweep (1- (array-dimension sweep 0)))
                haveSwept)))
              
          
(def resp  (+ (normal-rand (array-dimension poly 0))
              (matmult poly '(1 1 1 1))))


(print-matrix
 (covariance-matrix (bind-columns poly resp))  ; div by n-1
 t :float-digits 2)


(def sweep (make-sweep-matrix  poly resp))     ; X'X
(setf haveSwept ())
(print-matrix sweep t :float-digits 1) 


(def s (sweep 1))  (coefs)
(def s (sweep 2))  (coefs)
(def s (sweep 3))  (coefs)
(def s (sweep 4))  (coefs)



;--- Move the polynomials and see impact of collinearity

(setf poly (poly-matrix (+ 20 (rseq -1 1 100)) 4))
(def resp  (+ (normal-rand (array-dimension poly 0))
              (matmult poly '(1 1 1 1))))

(regression-model poly resp)

(def sweep (make-sweep-matrix  poly resp))     ; X'X
(setf haveSwept ())
(print-matrix sweep t :float-digits 2) 

(sweep 1)  (coefs)
(sweep 2)  (coefs)
(sweep 3)  (coefs)
(sweep 4)  (coefs)



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Decompositions: LU                                   ;;
;;                              Tierney, Section 5.5, 167  ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- LU factorization is most primitive, used by others (solve, detm, inv)
;    Result-> (LU in first matrix), pivoting, odd pivots, singular?

(lu-decomp '((1 2) (3 4)))       ; will not work with list

(setf lu (lu-decomp '#2a((1 2) (3 4))) )


;--- Solve once, use often

(lu-solve lu '(1 6))            

(lu-solve lu '(1 3))


;--- Works with complex numbers

(setf lu (lu-decomp '#2a((#c(1 0) #c(2 0))  (#c(1 1) #c(2 3)))  ))

(lu-solve lu '(1 3))




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Decompositions: Cholesky                             ;;
;;                              Tierney, Section 5.5, 168  ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Most often used with psd matrices. Result is (matrix val?) list

(def cov  (covariance-matrix
           (bind-columns (normal-rand 20) (normal-rand 20))))

(def CH (chol-decomp cov))

;--- Second element D is 0.0 if psd

ch

;--- Recover as M + D = L L'

(def c (first ch))

(print-matrix
 (- cov (matmult c (transpose c)))
 t :float-digits 2)


;--- Common application for correlated normal samples

   (defun CORR-NORMAL (n covMat)
     "Returns list of n normals with indicated covariance matrix."
     (let ((chol (chol-decomp covMat))  )
       (if (zerop (second chol))
           (let ((k (array-dimension covmat 0))  )
             (mapcar #'(lambda (x) (matmult (first chol) x))
                     (split-list (normal-rand (* n k)) k)))
           (format t "Cov not psd; offset by ~a~%" (second chol))
           )))

(def norm (corr-normal 50 '#2a((1 .5 0) (.5 1 .5) (0 .5 1) )))

(print-matrix 
 (setf cov (apply #'covariance-matrix (transpose norm)))
 t :float-digits 2)
     
       

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Decompositions: QR                                   ;;
;;                              Tierney, Section 5.5, 168  ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Embodies the Gram-Schmidt process (pivoting is optional)

(def data norm)
(setf qr (qr-decomp data))

(def data (make-array '(50 3) :initial-contents norm))
(setf qr (qr-decomp data))
          

;--- Note the form of the orthogonality ...

(print-matrix (covariance-matrix (first qr)) t :float-digits 2)

(def q (first qr))

(print-matrix (matmult (transpose q) q))


;--- Triangular structure (computed one row at a time)

(print-matrix (second qr))


;--- Relationship to Cholesky for square

(print-matrix (first (chol-decomp (matmult (transpose data) data))))

(print-matrix (second (qr-decomp data)))


              

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Decompositions: Eigenvalue                           ;;
;;                              Tierney, Section 5.5, 169  ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Get a list of evals, evectors

(print-matrix cov)

(def e (eigen cov))

(def L (first (chol-decomp cov)))


;--- Alternative square root

(print-matrix (- cov (matmult L (transpose L))))

 
    (defun MAT-SQRT (m)
      "Eigen square root of matrix"
      (let ((e (eigen m))  )
        (apply #'+
               (mapcar #'(lambda (val vec)
                           (* (sqrt val)
                              (outer-product vec vec #'*)))
                       (coerce (first e) 'list)
                       (second e)
                       ))
        ))


(def sr (mat-sqrt cov))
(print-matrix (- cov (matmult sr (transpose sr))))