;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Version 3/25/2000
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; methods for the arc proto to keep track of strange IR titles.
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defmeth arc :ir-title (&key title) title)
(defmeth arc :ir-short-title (&key short-title) short-title)
(defmeth arc :ir-method (&key method) method)
(defmeth arc :ir-response (&key response) response)
(defmeth arc :ir-titles () 
 (mapcar #'(lambda (a) (apply #'send arc :ir-title a)) (get-ir-meths)))
(defmeth arc :ir-short-titles () 
 (mapcar #'(lambda (a) (apply #'send arc :ir-short-title a)) (get-ir-meths)))
(defmeth arc :ir-methods () 
 (mapcar #'(lambda (a) (apply #'send arc :ir-method a)) (get-ir-meths)))
(defmeth arc :ir-responses () 
 (mapcar #'(lambda (a) (apply #'send arc :ir-response a)) (get-ir-meths)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;  Polynomial Inverse Regression
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Constant Variance Version of Polynomial inverse regression
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;(defparameter *ir-poly* nil)
;(defparameter *ir-nc-poly* nil)
(defparameter *ir-poly* t)
(defparameter *ir-nc-poly* t)
;; Add methods to dialog, unless they are already there


(defun get-ir-meths (&key (phdy *ir-phdy*) (poly *ir-poly*) (nc-poly *ir-nc-poly*))
"Function args: (&key (phdy *ir-phdy*) (poly *ir-poly*) (nc-poly *ir-nc-poly*))
Adds extra ir methods to the menu."
  (append *ir-meths*
   (if phdy
        (list (list :title "pHd(response)" :short-title "pHd-y"
              :method :phdy :response :yvar :test :phd-test)))
   (if poly
     (list (list :title "Quadratic fit" :short-title "IR2"
                 :method :ir-quad :response :yvar :test :poly-test)
           (list :title "Cubic fit" :short-title "IR3"
                 :method :ir-cubic :response :yvar :test :poly-test)
	   (list :title "Quartic fit" :short-title "IR4"
                 :method :ir-quartic :response :yvar :test :poly-test)
	   (list :title "Quintic fit" :short-title "IR5"
                 :method :ir-quintic :response :yvar :test :poly-test)))
   (if nc-poly
      (list (list :title "Quadratic-NC" :short-title "IRN2"
                 :method :ir-nc-quad :response :yvar :test :nc-poly-test)
           (list :title "Cubic-NC" :short-title "IRN3"
                 :method :ir-nc-cubic :response :yvar :test :nc-poly-test)
	   (list :title "Quartic-NC" :short-title "IRN4"
                 :method :ir-nc-quartic :response :yvar :test :nc-poly-test)
	   (list :title "Quintic-NC" :short-title "IRN5"
                 :method :ir-nc-quintic :response :yvar :test :nc-poly-test))))) 

;; Method to retrieve the degree of the polynomial
(defmeth inverse-regression-model-proto :degree () (slot-value 'degree))

;; Calling methods: set the degree and then do eigen computations in :ir-poly
(defmeth inverse-regression-model-proto :ir-quad () 
  (send self :add-slot 'degree 2) (send self :ir-poly))
(defmeth inverse-regression-model-proto :ir-cubic () 
  (send self :add-slot 'degree 3) (send self :ir-poly))
(defmeth inverse-regression-model-proto :ir-quartic () 
  (send self :add-slot 'degree 4) (send self :ir-poly))
(defmeth inverse-regression-model-proto :ir-quintic () 
  (send self :add-slot 'degree 5) (send self :ir-poly))

(defmeth inverse-regression-model-proto :ir-nc-quad () 
  (send self :add-slot 'degree 2) (send self :ir-nc-poly))
(defmeth inverse-regression-model-proto :ir-nc-cubic () 
  (send self :add-slot 'degree 3) (send self :ir-nc-poly))
(defmeth inverse-regression-model-proto :ir-nc-quartic () 
  (send self :add-slot 'degree 4) (send self :ir-nc-poly))
(defmeth inverse-regression-model-proto :ir-nc-quintic () 
  (send self :add-slot 'degree 5) (send self :ir-nc-poly))


;From eb@phalen.stat.umn.edu Wed Jun 26 14:39 CDT 1996
(defmeth inverse-regression-model-proto :sir-cov-mat-decomp ()
  (send self :covariance-matrix-decomp))


(defmeth inverse-regression-model-proto :pir-y (&key (centered nil) (all nil))
"Method args: (&key (centered nil) (all nil))
Returns the design matrix in PIR with deleted observations dropped, and centered if centered is t."
  (let* ((inc (which (send self :included)))
	 (k (send self :degree))
	 (y (select (send self :yvar) inc))
	 (design-mat (mapcar #'(lambda (pow) (^ y pow)) (iseq 1 k))))
    (if centered 
	     (- design-mat (mapcar #'mean design-mat))
             design-mat)))

(defmeth inverse-regression-model-proto :coefmat ()
"Returns the coefficient matrix of the p inverse regressions"

   (let* ((inc (which (send self :included)))
          (k (send self :degree))
          (x-centered (send self :sir-x :centered t))
          (x-std (matmult x-centered (send self :sir-cov-mat-decomp)))
          (x (transpose (array-to-nested-list x-std)))
	  (design-mat-ctr (apply #'bind-columns (send self :pir-y  :centered t)))
	  (x-mat (apply #'bind-columns x))
	  (pr (matmult (inverse (matmult (transpose design-mat-ctr) design-mat-ctr)) (transpose design-mat-ctr)))
	  (coefs (matmult pr x-mat))) coefs))

(defmeth inverse-regression-model-proto :z-mat ()
"Returns a square root of the left matrix multiplier of the 
coefficient matrix." 
 (let* ((k (send self :degree))
        (design-mat-ctr (apply #'bind-columns (send self :pir-y :centered t)))
	(n (send self :num-included))
        (d (matmult (transpose design-mat-ctr) design-mat-ctr))
        (sd (sv-decomp1 d)))
 (/ (matmult (first sd) (diagonal (^ (second sd) .5)) (transpose (third sd))) 
    (sqrt n))))


(defmeth inverse-regression-model-proto :parir-cov-mat-decomp ()
"Returns a square root of the inverse of the estimated covariance of X 
given Y based on the singular value decomposition." 
 (let* ((k (send self :degree))
        (n (which (send self :included)))
        (x-centered (send self :sir-x :centered t))
        (x-std (matmult x-centered (send self :sir-cov-mat-decomp)))
        (design-mat-ctr (apply #'bind-columns (send self :pir-y  :centered t)))
	(e (- x-std (matmult design-mat-ctr (send self :coefmat))))
        (m (send self :num-included))
        (sigma (/ (matmult (transpose e) e) m))
        (invsv (sv-decomp (inverse sigma))))
  (matmult (first invsv) (diagonal (^ (second invsv) .5)) 
            (transpose (third invsv)))))

(defmeth inverse-regression-model-proto :ir-poly ()
"Returns the eigenvalues of the standardized coefficient matrix,
and the eigenvectors of the space."  
 (let* ((k (send self :degree))
        (m (send self :num-included))
        (n (iseq (send self :num-cases)))
        (x (send self :x))
	(p (array-dimension x 1))
	(design-mat-ctr (apply #'bind-columns (send self :pir-y  :centered t)))
        (z-mat (send self :z-mat))
        (coefmat (send self :coefmat))
        (cov-mat (send self :parir-cov-mat-decomp))
        (b-stand (* (matmult z-mat coefmat cov-mat) (sqrt m)))
        (bsv (sv-decomp1 b-stand))
        (evalues (^ (second bsv) 2))
	(ord (reverse (order evalues)))
        (sp (matmult (transpose coefmat) 
                     (transpose design-mat-ctr) design-mat-ctr coefmat))
        (evectors (select (eigenvectors sp) ord))
        (invsigmax (sv-decomp (inverse (covariance-matrix x))))
        (invhalf (matmult (first invsigmax) 
           (diagonal (^ (second invsigmax) .5)) (transpose (third invsigmax))))
	(evect-to-list (mapcar #'(lambda (v) (coerce (matmult invhalf v) 'list)) evectors))
)   
    (setf (slot-value 'eigen-structure) 
    (list (coerce (select evalues ord) 'list)
       (select evect-to-list ord)))))

(defun sv-decomp1 (x)
"Function args: (x)
Like sv-decomp, except it allows the n by p matrix x to have n < p."
 (let* ((dim (array-dimensions x)))
 (cond ((apply #'< dim)
        (let ((ans (sv-decomp (transpose x))))
         (list (transpose (third ans)) (second ans)
               (transpose (first ans)) (select ans 3))))
       (t (sv-decomp x)))))

(defmeth inverse-regression-model-proto :poly-test ()
"Message args: ()
Returns the p partial sums of ranked eigenvalues of the standardized 
coefficient matrix that can be compared to the
percentage points of a chi-square distribution with (q-s)(p-s) degrees of
freedom, where s = number of terms in the partial sum, 
s= 1,2,...,min(q,p), p=number of x variables, q=k, where 
k=polynomial degree."
  (let* ((k (send self :degree))
         (eigen (slot-value 'eigen-structure))
         (e (reverse (first eigen)))
	 ;(eigen (select (send self :ir-poly) 0))
	 ;(e (reverse eigen))
         (p (array-dimension (send self :x) 1))
         (q k)
         (nobs (send self :num-included)))
  (when (> k 0)
   (format t "~%Approximate Chi-squared test statistics based on partial sums")
   (format t "~%of singular values times ~a~%~%" nobs)
   (format t "Number of   ~13tTest~%")
   (format t "Components  ~13tStatistic~27tdf~30tp-value~%")
   (dolist (i (reverse (iseq 0 (- (min p q) 1))))
         (let ((a (sum (select e (iseq 0 i))))
               (b (- 1 (chisq-cdf (sum (select e (iseq 0 i)))
                              (* (- q (- (- (min q p) 1) i)) 
                                 (- p (- (- (min q p) 1) i))))))
               (c (* (- q (- (- (min q p) 1) i)) (- p (- (- (min q p) 1) i)))))
        (format t "~4d ~12t~11a~22t~5d~32t~5,3f~%" (- p i) (optfmt1 a) c b))))))




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Non - Constant Variance Version of Polynomial inverse regression
;;; 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; The following w-mat method returns the Wn=(Fn'Fn)^-1Fn' matrix
;; in "Estimating the Structural Dimension of Regressions via
;; parametric inverse regression"

;; The following w-mat method returns the Wn=(Fn'Fn)^-1Fn' matrix
;; in the paper above

(defmeth inverse-regression-model-proto :w-mat ()
"Returns the left matrix multiplier of the regressor matrix in the formula of
the coefficient matrix." 
 (let* ((k (send self :degree))
        ;(y (select (send self :yvar) (which (send self :included))))
        ;(design-mat (apply #'bind-columns 
        ;   (mapcar #'(lambda (pow) (^ y pow)) (iseq (1+ k)))))
        (design-mat-ctr (apply #'bind-columns (send self :pir-y :centered t)))
	(n (send self :num-included))
        (d (matmult (transpose design-mat-ctr) design-mat-ctr))
        (invd (inverse d)))
    (matmult invd (transpose design-mat-ctr))))


(defmeth inverse-regression-model-proto :nc-cov-mat ()
"Returns the pq x pq estimated covariance matrix in paper2. That is the ijth block is given by Wn diagonal (sigma_ij(Y)) Wn' (the diagonal is of dimension n)"
 (let*  ((k (send self :degree))
         (p (array-dimension (send self :x) 1))
         (q k)
         (pq (* p q))
	 (inc (which (send self :included)))
         (w (send self :w-mat))
	 (x-centered (send self :sir-x :centered t))
	 (x-std (matmult x-centered (send self :sir-cov-mat-decomp)))
	 (z (transpose (array-to-nested-list x-std)))
	 (y (select (send self :yvar) inc))
         (H (make-array (list pq pq) :element-type 'c-double)))
    (flet
      ((sblock (i j) (matmult w (diagonal (- (send (regression-model y (* (select z i) (select z j)) :print nil) :fit-values) (* (send (regression-model y (select z i) :print nil) :fit-values) (send (regression-model y (select z i) :print nil) :fit-values)))) (transpose w))))
   
      ;; Form H block-wise
      ;; H is a p-by-p array of q-by-q sym. matrices
      ;; Use symmetry to avoid unnecessary work
      (do ((i 0 (+ i 1))
           (r (iseq q) (+ r q)))
          ((= i p) H)
          
          ;; diagonal block
          (setf (select H r r) (sblock i i))
          ;; off-diagonal blocks
          ;; below and right of diagonal block
          (do ((j (+ i 1) (+ j 1))
               (c (+ r q) (+ c q)))
              ((= j p) nil)
              (setf (select H r c) (sblock i j))
              (setf (select H c r) (select H r c)))))))



(defmeth inverse-regression-model-proto :ir-nc-cov-mat-decomp ()
"Returns a square root of the inverse of the estimated non-constant covariance  of X given Y based on the singular value decomposition. See paper2 for
 details." 
 (let* ((k (send self :degree))
	(m (send self :num-included))
	(sigma (* m (send self :nc-cov-mat)))
	;(sigma (send self :nc-cov-mat ))
	(invsv (sv-decomp (inverse sigma))))
  (matmult (first invsv) (diagonal (^ (second invsv) .5)) 
            (transpose (third invsv)))))

(defmeth inverse-regression-model-proto :ir-nc-poly ()
"Returns the eigenvalues of the standardized coefficient matrix,
and the eigenvectors of the space."  
 (let* ((k (send self :degree))
        (m (send self :num-included))
        (n (iseq (send self :num-cases)))
	(q k)
        ;(x (send self :sir-x :centered t))
        (x (send self :x))
	(p (array-dimension x 1))
	;(y (select (send self :yvar) (which (send self :included))))
        ;(design-mat (apply #'bind-columns 
        ;    (mapcar #'(lambda(pow) (^ y pow)) (iseq (1+ k)))))
	(design-mat-ctr (apply #'bind-columns (send self :pir-y  :centered t)))
	(h-mat (send self :ir-nc-cov-mat-decomp))
        (coefmat (send self :coefmat))
	(coefmat-to-vector (make-array (* p q) :displaced-to coefmat))
        ;(b-stand-vec (* (matmult h-mat coefmat-to-vector) (sqrt m)))
	(b-stand-vec (matmult h-mat coefmat-to-vector))
	(dim (list q p))
	(b-stand (make-array dim :displaced-to b-stand-vec))
        (bsv (sv-decomp1 b-stand))
        ;(evalues (^ (second bsv) 2))
	(evalues (* (^ (second bsv) 2) m))
	;(ord-eval (iseq (min p q)))
        (ord (reverse (order evalues)))
	(sp (matmult (transpose coefmat) 
                     (transpose design-mat-ctr) design-mat-ctr coefmat))
        (evectors (select (eigenvectors sp) ord))
        (invsigmax (sv-decomp (inverse (covariance-matrix x))))
        (invhalf (matmult (first invsigmax) 
           (diagonal (^ (second invsigmax) .5)) (transpose (third invsigmax))))
	(evect-to-list (mapcar #'(lambda (v) (coerce (matmult invhalf v) 'list)) evectors))
)   (setf (slot-value 'eigen-structure) 
    (list (coerce (select evalues ord) 'list)
       (select evect-to-list ord)))))


(defmeth inverse-regression-model-proto :nc-poly-test ()
"Message args: ()
Returns the p partial sums of ranked eigenvalues of the standardized 
coefficient matrix that can be compared to the
percentage points of a chi-square distribution with (q-s)(p-s) degrees of
freedom, where s = number of terms in the partial sum, 
s= 1,2,...,min(q,p), p=number of x variables, q=k, where 
k=polynomial degree."
  (let* ((k (send self :degree))
         (eigen (slot-value 'eigen-structure))
         (e (reverse (first eigen)))
	 (p (array-dimension (send self :x) 1))
         (q k)
         (nobs (send self :num-included)))
  (when (> k 0)
   (format t "~%Approximate Chi-squared test statistics based on partial sums")
   (format t "~%of singular values times ~a~%~%" nobs)
   (format t "Number of   ~13tTest~%")
   (format t "Components  ~13tStatistic~27tdf~30tp-value~%")
   (dolist (i (reverse (iseq 0 (- (min p q) 1))))
         (let ((a (sum (select e (iseq 0 i))))
               (b (- 1 (chisq-cdf (sum (select e (iseq 0 i)))
                              (* (- q (- (- (min q p) 1) i)) 
                                 (- p (- (- (min q p) 1) i))))))
               (c (* (- q (- (- (min q p) 1) i)) (- p (- (- (min q p) 1) i)))))
        (format t "~4d ~12t~11a~22t~5d~32t~5,3f~%" (- p i) (optfmt1 a) c b))))))