#|

     class12.lsp   

    - Better moving average algorithms
    
    - Timing comparisons:  is the time for foreign functions well-spent?


|#


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Moving average simulation                            ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;--- 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 p (plot-points xax (setf DATA (sine-gen n noiseSD)) ))
(send p :add-lines xax truth :color 'magenta)


;--- Moving average estimator

   (defun MOVING-AVERAGE-1 (x width)
     "Smooths assuming input x is equally spaced."
     (unless (oddp width)
             (format t "Width needs to be odd integer.~%")
             (setf width (let ((rw (round width))  )
                           (if (oddp rw) rw (1+ rw))))   )
     (let ((half (/ (1- width) 2))      ; invariants done once
           (n-1  (1- (length x)))
           (smth ())  )
       (flet ((avg (i) 
                   (mean (select x
                                 (iseq (max 0   (- i half))
                                    (min n-1 (+ i half)))))  ))
         (dotimes (i (length x))
                  (push (avg i) smth))
         (nreverse smth)
         )))

(moving-average-1 '(1 2 3 4 5) 3)


;--- 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)))


(moving-average '(1 2 3 4 5) 3)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Timing comparisons                                   ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(time
 (dotimes (j 100)
          (moving-average-1 (iseq 100) 3)))

The evaluation took 5.03 seconds; 0.10 seconds in gc. ; n=100, len  3
The evaluation took 6.87 seconds; 0.40 seconds in gc. ;      , len 15
The evaluation took 8.52 seconds; 0.45 seconds in gc. ;      , len 31


;--- now compiled

(compile 'moving-average-1)

The evaluation took 4.12 seconds; 0.45 seconds in gc. ;      , len 31


;--- now better alg, no respose to window length

(time
 (dotimes (j 100)
          (moving-average   (iseq 100) 31)))

The evaluation took 5.50 seconds; 0.17 seconds in gc. ; n=100, len  3
The evaluation took 5.33 seconds; 0.08 seconds in gc. ;            15
The evaluation took 5.27 seconds; 0.13 seconds in gc. ;            31
 
(compile 'moving-average)

The evaluation took 0.77 seconds; 0.05 seconds in gc.


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


;--- Fun with animation, slider

(def sin (sample (iseq n) (min n 100)))         ; show at most 100 pts
(def p2 (plot-points (select xax sin)
                     (select data sin)))  


(interval-slider-dialog            
 '(10 500)       ; limits on interval
 :points 50      ; divisions of interval
 :action (lambda (x)
           (send p2 :clear-lines)
           (send p2 :add-lines xax (moving-average data x))
           ))
         


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;     Foreign functions                                   ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Compile the C code, then prep it as a shared library

cc -c +O3 -Aa +Z  movAvg.c

ld -o movAvg.shlb -b movAvg.o


;--- Load the C function library into XLisp

(dyn-load "/users/bob/html/stat540/movAvg.shlb")


;(open-resource-file "Dave:Text:Classes:Stat540:movAvg.shlb"
;                    :verbose t)

;--- call the C function (case matters!)

(call-cfun name arg1 arg2 ... )

    (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))))   )
      (call-cfun "MovAvg"
                 (float x) (length x)
                 width
                 (repeat 7.7 (length x))
                 )
      )

             
(mov-avg '(1 2 3 4 5) 3)


(time
 (dotimes (j 100)
          (mov-avg   (iseq 100) 31)))
The evaluation took 0.19 seconds; 0.01 seconds in gc.  ; worth it???


; --- difference emerges as x gets longer...

(def x (uniform-rand 1000))

(time (dotimes (i 10)
               (moving-average x 151)))   ; 3.48 sec

(time (dotimes (i 10)
               (mov-avg x 151)))           ; 0.68 Mac, 0.17 HP



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;     Basic array manipulation                            ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Create an array in various ways:  key, explicit, from lists, indirectly

(def m '#2a((1 2) (3 4)))

(def m (make-array '(3 4) :initial-contents (iseq 12)))    ; initial-element

(def m (bind-rows '(1 2 3) '(4 5 6)))                      ; bind columns

(def iden (identity-matrix 5))

(def diag (diagonal '(4 5 6)))

(def cov (covariance-matrix
	 (apply #'bind-columns
	 	(mapcar #'(lambda (n) (normal-rand n))  '(100 100 100)))))
	 	
 	
;--- Conversions

(coerce '(1 2 3 4) 'vector)

(coerce '#2a((1 2 3 4 5)) 'vector)

(compound-data-seq '#2a((1 2 3) (4 5 6)))                  ; row major order



;--- Display

(print-matrix cov)

(print-matrix cov t :float-digits 4)

	 	
;--- Dimensions, "shape" of an array

(array-dimensions m)

(array-dimension m 0)

(rank m)  ; bugged in this version

(transpose m)


;--- Indexing

(aref cov 1 1)

(select cov nil 1)

(select cov (iseq 3) 1)


;--- Operations

(matmult diag cov)

(print-matrix 
	(matmult cov (inverse cov))  )


(solve cov '(1 0 0))

(inverse cov)

(determinant cov)

(setf ee (eigen cov))

(matmult cov (/ (first (second ee)) (aref (first ee) 0)))


;--- Sweep matrix

(def x1 (normal-rand 10))
(def x2 (normal-rand 10))

(def y  (+ x1 (* -1 x2) (* .5 (normal-rand 10))))

(def rm (regression-model (list x1 x2) y))


(def sweep 
	(make-sweep-matrix (bind-columns x1 x2) y))  ; opt weights

(print-matrix sweep t :float-digits 2)

(def swept
	(sweep-operator sweep '(1 2)))

(print-matrix (first swept) t :float-digits 2)


;--- See Tierney, Section 5.5 for more hints on how to use the Sweep function