;;;
;;;  General Partial SIR version of March 25, 2005
;;; 
;;;
;;; Menu items
;;;

(when (null (position 'fit-gpsir-menu-item *arc-fit-menu-items* :test #'eq))
 (defparameter *arc-fit-menu-items*
   (append *arc-fit-menu-items* (list 'fit-gpsir-menu-item))))
(rc-menu-item 'fit-gpsir-menu-item     "Fit GPartial SIR ...      " :fit-gpartial-sir)


;;;
;;; Define the prototype
;;;

(defproto gpartial-sir-proto
  '(y1 response1-name response2-name y2 nslice1 nslice2 
    actual-nslices1 actual-nslices2 gpdim gpdir) 
  ()
  inverse-regression-model-proto)

(defmeth gpartial-sir-proto :isnew (&rest args &key y1 response1-name y2 
    response2-name nslice1 nslice2 actual-nslices1 actual-nslices2 gpdim gpdir)
  
  (send self :y1 y1)
  (send self :y2 y2)  
  (send self :method :gpsir) 
  (send self :method-response :y)
  (send self :response1-name response1-name)
  (send self :response2-name response2-name)
  (send self :nslice1 nslice1)
  (send self :nslice2 nslice2)
  (send self :gpdim gpdim)
  (apply #'call-next-method args))

(defmeth gpartial-sir-proto :nslices (&rest args) (+ (* (send self :nslice1) 
(send self :nslice2)) 1))

(defmeth gpartial-sir-proto :num-cases () (length (send self :y1)))


;;;
;;; In the next method :y the method response is Y2
;;;

(defmeth gpartial-sir-proto :y (&rest args) 
(send self :y2))

(defmeth gpartial-sir-proto :y1 (&optional new-y1)
"Method args: (&optional new-y1)
With no argument returns the y sequence as supplied to m. With an argument
NEW-Y sets the y sequence to NEW-Y and recomputes the estimates. Y1 plays 
the role of the conditioning variable W in the paper."
  (when (and new-y1 (or (matrixp new-y1) (sequencep new-y1)))
        (setf (slot-value 'y1) new-y1)
        (send self :needs-computing t))
  (slot-value 'y1))


(defmeth gpartial-sir-proto :y2 (&optional new-y2)
"Method args: (&optional new-y2)
With no argument returns the y sequence as supplied to m. With an argument
NEW-Y sets the y sequence to NEW-Y and recomputes the estimates. This is 
the response Y in the paper."
  (when (and new-y2 (or (matrixp new-y2) (sequencep new-y2)))
        (setf (slot-value 'y2) new-y2)
        (send self :needs-computing t))
  (slot-value 'y2))

(defmeth gpartial-sir-proto :response1-name (&optional name1)
  (if name1 (setf (slot-value 'response1-name) name1))
  (slot-value 'response1-name))

(defmeth gpartial-sir-proto :response2-name (&optional name2)
  (if name2 (setf (slot-value 'response2-name) name2))
  (slot-value 'response2-name))

(defmeth gpartial-sir-proto :nslice1 (&optional new-n)
"Method args: (&optional new-n)
With no argument returns the number of slices used by pSIR for response 1,
that is W.  With an argument NEW-N sets the number of slices 
to NEW-N and recomputes the estimates"
  (when (and new-n (> new-n 0))
        (setf (slot-value 'nslice1) new-n)
        (send self :needs-computing t))
  (slot-value 'nslice1))

(defmeth gpartial-sir-proto :nslice2 (&optional new-n) 
"Method  args: (&optional new-n) With no argument returns the number of 
slices used by pSIR for response 2, that is Y.  With an argument 
NEW-N sets the number of slices to NEW-N and recomputes the 
estimates"
  (when (and new-n (> new-n 0))
        (setf (slot-value 'nslice2) new-n)
        (send self :needs-computing t))
  (slot-value 'nslice2))

(defmeth gpartial-sir-proto :gpdim (&optional new-n)
"Method args: (&optional new-n)
With no argument returns the dim d used by gpSIR.  
With an argument NEW-N sets the d to NEW-N and recomputes the
estimates"
  (when (and new-n (> new-n 0))
        (setf (slot-value 'gpdim) new-n)
        (send self :needs-computing t))
  (slot-value 'gpdim))
  
(defmeth gpartial-sir-proto :gpdir ()
"Method args: ()
With no argument returns the d-directions produced by gpSIR."
    (if (send self :needs-computing) (send self :compute))
  (slot-value 'gpdir))
  
;;;
;;; GPartial sir dialog
;;;

(defun gpartial-sir-dialog (names &key (name "GPartial SIR") (weights "")
          (response1 "") nslice1 (response2 "") nslice2  gpdim
          (candidates names) (predictors (repeat "" (length names))))
"Function args: (names &key name)"
  (let* ((version (send text-item-proto :new
         (format nil "R-code ~a" *arc-version*)))
         (all (send allocate-items-proto :new))
         (candidates-l (send subordinate-list-proto :new 
                        (pad-list candidates (length names)) 'primary all))
    (candidates-n (send text-item-proto :new "Candidates"))
    (predictors-l (send subordinate-list-proto :new 
                        (pad-list predictors (length names)) 'secondary all))
    (predictors-n (send text-item-proto :new "Predictors"))
    (response1-l (send subordinate-list-proto :new (list response1)
                        'oneitem all))
    (response1-n (send text-item-proto :new "Condition on W...  "))
    (response2-l (send subordinate-list-proto :new (list response2)
                        'oneitem all))
    (response2-n (send text-item-proto :new "Response Y...  "))
    (weight-l (send subordinate-list-proto :new 
          (if weights (list weights) '("")) 'oneitem all))
    (weight-n (send text-item-proto :new "Weights...     "))
    (mod-type (format nil "GPartial SIR"))
    (name-l (send edit-text-item-proto :new name :text-length 14 ))
    (name-n (send text-item-proto :new 
             (format nil "Name for ~%~a" mod-type)))
    (slice-n (send text-item-proto :new
             (format nil "Number of~%Slices for")))
    (nslice1-l (send edit-text-item-proto :new
             (format nil "~a" nslice1) :text-length 4))
    (nslice1-n (send text-item-proto :new
             (format nil "W")))
    (nslice2-l (send edit-text-item-proto :new
             (format nil "~a" nslice2) :text-length 4))
    (nslice2-n (send text-item-proto :new (format nil "Y")))
    (gpdim-l (send edit-text-item-proto :new
             (format nil "~a" gpdim) :text-length 4))
    (gpdim-n (send text-item-proto :new
             (format nil "d")))
    (get-answer
     #'(lambda ()
         (let* ((ans (send all :value))
                (raw (send all :all-values))
                (slc1 (with-input-from-string (s (send nslice1-l :text))
                                              (read s nil)))
                (slc2  (with-input-from-string (s (send nslice2-l :text))
                                              (read s nil))) 
                (gpd (with-input-from-string (s (send gpdim-l :text))
                                              (read s nil))))                                             
                                              ;end lets
           (cond ((or (null (select ans 1))
                      (< (length (select ans 1)) 2)
                      (null (select ans 2))
                      (null (select ans 3))
                      (not (integerp slc1))
                      (not (integerp slc2))
                      (not (integerp gpd))
                      (not (plusp slc1))
                      (not (plusp slc2))
                      (not (plusp gpd)))
                  (message-dialog 
                    "Predictors or responses missing~%or badslice numbers")
                  (gpartial-sir-dialog names
                            :name (send name-l :text)
                            :candidates (select raw 0)
                            :predictors (select raw 1)
                            :response1 (first (select raw 2))
                            :response2 (first (select raw 3))
                            :weights (first (select raw 4))
                            :nslice1 slc1
                            :nslice2 slc2
                            :gpdim gpd))
                  (t (list :predictors (select ans 1)
                            :candidates (select ans 0)
                            :response1 (first (select ans 2))
                            :response2 (first (select ans 3))
                            :weights (first (select ans 4))
                            :name (send name-l :text)
                            :nslice1 slc1
                            :nslice2 slc2
                            :gpdim gpd )))))) ;end get-answer
    (cancel (send modal-button-proto :new "Cancel"))
    (done (send modal-button-proto :new "Done" :action 
                     #'(lambda () (funcall get-answer))))
    (dialog (send modal-dialog-proto :new
                       (list (list version name-n name-l)
                             (list (list candidates-n candidates-l 
                response1-n response2-n
                weight-n)
               (list predictors-n predictors-l 
                response1-l response2-l
                weight-l)
                                   (list slice-n
                (list nslice1-n nslice1-l)
                (list nslice2-n nslice2-l)
                (list gpdim-n gpdim-l)
                done cancel))))))
   (send dialog :default-button done)
   (send dialog :modal-dialog) ))

;;;
;;; Determine slices
;;;

(defun gpartial-sir-slices (y1 y2 num-slices1 num-slices2)
  (let* ((first-slice (gpartial-sir-first-slice y1 num-slices1)))
    (apply #'append (gpartial-sir-second-slice y2 first-slice num-slices2))))

(defun gpartial-sir-first-slice (y num-slices)
  (let  ((ind y))
    (sir-slices ind num-slices)))

(defun gpartial-sir-second-slice (y2 first-slice num-slices)
  (mapcar #'(lambda (fs)
         (let  ((y2i (select y2 fs)))
      (mapcar #'(lambda (inds)
             (select fs inds))
         (sir-slices y2i num-slices)) ))
     first-slice))
     
;;;;
;;; GPartial SIR methods
;;;;

(defmeth gpartial-sir-proto :class-slice-indices (n)
"Message args: (n)
Returns as a list of lists the indices relative to the included data over 
all classes of the slices of Y2 corresponding to the slices within slice n of Y1.  
If cases have been deleted then indices are relative to the reduced data set. Used for 
general partial SIR when Y1 is the conditioning predictor."
  (let* ( 
         (inc (which (send self :included)))
         (y11 (select (send self :y1) inc))
         (y21 (select (send self :y2) inc))
         (num-slices1 (+ (send self :nslice1) 1)) 
; 1 was added in the last statement to insure proper no.
         (num-slices2 (send self :nslice2))
         (first-slice (gpartial-sir-first-slice y11 num-slices1)))
    (first (gpartial-sir-second-slice y21 (list (nth n first-slice)) num-slices2))))

(defmeth gpartial-sir-proto :class-num-obs ()
"Message args: ()
Returns a list of the number of observations in the slices of Y1 (W), the 
conditioning variable."
(let* ( 
        (inc (which (send self :included)))
        (y11 (select (send self :y1) inc))
        (num-slices1 (+ (send self :nslice1) 1)) 
; 1 was added in the last statement to ensure proper no.
        (first-slice (gpartial-sir-first-slice y11 num-slices1))
        (wt (mapcar #'length first-slice))
)
wt))

(defmeth gpartial-sir-proto :class-wts ()
"Message args: ()
Returns a list of the weights n_c/n for combining over the classes c of 
the conditioning variable Y1 (W). n is 
the total no of obs and n_c is the no of obs in class c of W"
(let* ( 
       (nc (send self :class-num-obs))
(n (sum nc))
)
(/ nc n)))

(defmeth gpartial-sir-proto :class-y2-slices (ic)
"Message args: (n)
Returns a list of the indices of the sliced response y2
 (not conditioning variable y1) in class ic.  The indices are relative to just the 
included data in class ic."
(let* (
       (inc (which (send self :included)))
       (y11 (select (send self :y1) inc))
       (y21 (select (send self :y2) inc))
       (nslc1 (+ (send self :nslice1) 1))
       (class-ind (nth ic (gpartial-sir-first-slice y11 nslc1)))
       (y211 (select y21 class-ind))
       )
  (sir-slices y211 (send self :nslice2))
  ))

(defmeth gpartial-sir-proto :class-x (ic &key (centered nil))
"Message args: (ic)
Returns a list of 2 items: (1) the x matrix of the predictors corresponding to class n; 
that is, slice n of Y1, and centered if centered is t. (2) The corresponding case 
weights ."
(let* ((p (send self :basis))
       (inc (which (send self :included)))
       (wts (send self :sir-wts)) ;accounts for deleted cases.
       (wts (if wts wts (repeat 1 (length inc))))
       (y11 (select (send self :y1) inc))
       (nslc1 (+ (send self :nslice1) 1))
       (x (select (send self :x) inc p))
       (class-ind (nth ic (gpartial-sir-first-slice y11 nslc1)))
       (x1 (select x class-ind p))
       (wts (select wts class-ind))
       (means (mapcar #'mean (column-list x1)))
       (x1 (if centered 
           (apply #'bind-columns (mapcar #'- (column-list x1) means))
           x1))
       )
(list x1 wts)))

(defmeth gpartial-sir-proto :class-slice-means (ic)
"Method args: (c)
Returns the noslices times p matrix of slice means in the X scale 
(x-bar-i minus x-bar) for class c. The slice means are the rows of the 
matrix."
  (let* ((col-means #'(lambda (m) (sweep m 'cols #'(lambda (c) (mean c)))))
         (x (nth 0 (send self :class-x ic :centered t)))
         (slices (send self :class-y2-slices ic))
         (all (iseq (array-dimension x 1))))
    (apply #'bind-rows
             (mapcar #'(lambda (z) (funcall col-means (select x z all)))
                       slices))))

(defmeth gpartial-sir-proto :class-cov-means (c)
"Method args: (c)
Returns the covariance matrix of the slice means in the X scale for class c."
  (let* (
         (w (diagonal (mapcar #'length (send self :class-y2-slices c))))
         (means (send self :class-slice-means c))
         (z-trans (transpose means))
         (cov (/ (matmult (matmult z-trans w) means) (sum w)))
         )
    cov))

;;;
;;; Determine actual slices
;;;

(defmeth gpartial-sir-proto :actual-nslices1 ()
"Message args: ()
This returns the actual number of first slices."
  (setf (slot-value 'actual-nslices1) 
        (length (gpartial-sir-first-slice (send self :y1) (send self :nslice1))))
  (slot-value 'actual-nslices1))

(defmeth gpartial-sir-proto :actual-nslices2 ()
"Message args: ()
This returns the actual number of first slices."
  (setf (slot-value 'actual-nslices2) 
        (mapcar (function length) 
                (gpartial-sir-second-slice 
                 (send self :y2) 
                 (gpartial-sir-first-slice (send self :y1) (send self :nslice1))
                                     (send self :nslice2))))
  (slot-value 'actual-nslices2))

(defmeth gpartial-sir-proto :actual-nslices ()
(sum (send self :actual-nslices2)))
      
;;;;
;;;; unique for gp.sir
;;;;

(defmeth gpartial-sir-proto :gpsir ()
"Method args: ()
Puts a list of the directions in x-scale into the slot 
eigenstructure."
  (let* ((dim (send self :gpdim))
         (egn (send self :gpsir-alsm dim)))
    (setf (slot-value 'slices) (send self :sir-slices))
    (setf (slot-value 'gpdir) (second egn))
   ))

  
(defmeth gpartial-sir-proto :class-cov-matrix (c)
"Method args: (c)
Returns a list with the first element 
the covariance matrix for class c in the original scale.  The second is 
the number of observations making up that matrix.  Weights 
and deleted cases OK"
(let* (
(x (first (send self :class-x c)))
(cov (covariance-matrix x))
)
(list cov (first (array-dimensions x)))
))


(defmeth gpartial-sir-proto :test () (send self :gpsir-test))

(defmeth gpartial-sir-proto :gpsir-test 
    (&key (dim 3) (print t) )
"Method args: (&key (print t) )
Returns the MDF values which can be compared the 
percentage points of a chi-square distribution with a linear 
combination of iid chi-square one variates."
(let* (  (h (send self :actual-nslices))
         (c (send self :actual-nslices1))
         (p (length (send self :basis)))
         (hmp (- h c p)) 
         (nobs (send self :num-included))
        a temp coefs  avep (pvals ()) (flag ()))
  (when print
    (format t "~%Approximate Chi-squared test statistics")
    (format t "~%Number of   ~13tTest~%")
    (format t "Components  ~13tStatistic~24ttrace~32tp-value~28tflag~%")
    (do ((i 0 (+ i 1))) ((= i (min (- h 1) p dim)))
      (let* ((temp (send self :gpsir-alsm i))        
             (a (nth 0 temp))
             (coefs (send self :gpsir-coefs i (nth 1 temp) (nth 2 temp))))                      
            (when (> a 0)
              (setf avep (/ (+ (wood-pvalue (nth 0 coefs) a)
                 (- 1 (imhof-cdf a (nth 0 coefs)))) 2))
;;;              (format t "~4d ~12t~11a~22t~5d~32t~5,3f~40t~3d~%"
;;;                        (+ 1 i) a (sum (nth 0 coefs))
;;;                    (wood-pvalue (nth 0 coefs) a) (nth 1 coefs))
              (format t "~4d ~12t~11a~22t~5d~32t~5,3f~40t~3d~%"
                        (+ 1 i) a (sum (nth 0 coefs)) avep (nth 1 coefs))
              ))))
  (if (null print)
      (do ((i 0 (+ i 1))) ((= i (min (- h 1) p dim)))
        (let* ((temp (send self :gpsir-alsm i))        
           (a (nth 0 temp))
           (coefs (send self :gpsir-coefs i (nth 1 temp) (nth 2 temp))))
      (when (> a 0)
        (setf avep (/ (+ (wood-pvalue (nth 0 coefs) a) (- 1 (imhof-cdf a (nth 0 coefs)))) 2))
        (setf pvals (append pvals (list avep)))
;;;        (setf pvals (append pvals (list (wood-pvalue (nth 0 coefs) a))))
;;;       (setf pvals (append pvals (list (- 1 (imhof-cdf a (nth 0 coefs))))))    
            ))))
  (if (null print) pvals)
  )
)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;  gpsir-gammaw      

(defmeth gpartial-sir-proto :gpsir-gammaw (c)
"Method args: (c)
Return the covariance matrix gammaw of the c-th population of Sigmaw^{-1}\xbar_{wy}.
Notice it is always for d>0."
 (let* (
        (nw (nth c (send self :class-num-obs)))
        (sigmawinv (inverse  (first (send self :class-cov-matrix c))))
        (mxw (send self :class-slice-means c))
        (xw (first (send self :class-x c :centered T)))
        (slicesw (send self :class-y2-slices c))
        (hw (length slicesw))
        (p (length (send self :basis)))
        (nkw (mapcar 'length slicesw))
    (w (send self :sir-wts))
    (sw (mapcar #' (lambda (s) (sum (select w s))) slicesw))
        (gammaw 0) xjw veck contr slicekw temp )
;;; begin outer loop. 
    (dotimes (k hw) 
      (setf veck (repeat 0 hw))
      (setf (nth k veck) 1)
      (setf veck (- veck (/ nkw nw)))
      (setf slicekw (nth k slicesw))    
;;; begin inner loop. 
      (dotimes (j (nth k nkw))
;;;; here xjw is sigma^{-1}x.
    (setf xjw (matmult (select xw (nth j slicekw) 
                   (iseq p)) sigmawinv))
        (setf contr (matmult (outer-product xjw xjw) 
              (transpose mxw) (diagonal (/ sw nw))))
    (setf contr (- (outer-product xjw veck) contr))
;;; notice the transformation for :gpsir-coefs. 
    (setf contr (matmult  contr (diagonal (/ nw sw))))
    (setf contr (vecc contr))
        (setf gammaw (+ gammaw (outer-product contr contr)))
       ))
;;; end of inner and outer loop. 
    (/ gammaw nw)
))

;;;;;;;;;;;;;;;;;;;;;
;;;  gpsir-coefs  

(defmeth  gpartial-sir-proto :gpsir-coefs (d b0 c0)
 "Method args: (num-slices d b0 c0)
 Return the coefficients for the linear combination 
  of chi-squares when PCS is d-dimensional."
(let* (  (h (send self :actual-nslices))
         (numclass (send self :actual-nslices1))
         (p (length (send self :basis)))
         (vecnw (send self :class-num-obs))
     coefs cw sigmawhalf slicesw hw nkw fw 
     ctransp lmat  delta  u
     (middle ())  gammahalf (vhalf ())
     (rankdelta (* d (+ h p (- d))))      
     temp (flag 0))

  (if (= d 0) 
    (setf coefs (repeat 1 (* p (- h numclass))))
   (tagbody 
    (dotimes (w numclass)
        (setf sigmawhalf (sqroot-matrix (nth 0 (send self :class-cov-matrix w))))
        (setf slicesw (send self :class-y2-slices w))
        (setf hw (length slicesw))
        (setf nkw (mapcar 'length slicesw))
        (setf fw (/ nkw (sum nkw)))
        (setf cw (matmult (select c0 w) (diagonal (/ (sqrt nkw)))))
        (if (eq w 0) (setf ctransp cw)
            (setf ctransp (bind-columns ctransp cw)))
        (setf vhalf (diagonal-block vhalf 
        (kronecker-product (diagonal (sqrt fw)) sigmawhalf)))
        (setf gammahalf (diagonal-block gammahalf 
               (sqroot-matrix (send self :gpsir-gammaw w)))))    
    (setf delta (bind-columns (kronecker-product 
    (transpose ctransp) (identity-matrix p)) 
            (kronecker-product (identity-matrix h) b0)))
    (setf temp (eigen (matmult vhalf delta (transpose delta) vhalf)))
    (if  (nth 2 temp) (setf flag 1))
    (setf lmat (select (nth 1 temp) (iseq rankdelta)))
    (setf lmat (transpose (make-array (list rankdelta (* p h))
                :initial-contents (combine lmat))))
    (setf u (matmult vhalf (- (identity-matrix (* p h))
               (matmult lmat (transpose lmat))) vhalf))
    (setf coefs (eigenvalues  (matmult gammahalf u gammahalf))))
    )

 (list (coerce (select coefs (which (< 0.01 coefs))) 'list) flag)
))

 ;;;;;;;;;;;;;;;;;;;;
 ;;; gpsir-alsm   ;;;
 
(defmeth gpartial-sir-proto :gpsir-alsm (d)
"Method args: (d)
Return d directions which maximize the projection function
 and function value via iterative least square methods."
(let* (
        (itmax 1000) (eps 1.0E-6)
        (ni (send self :class-num-obs))
        (p (length (send self :basis)))
        (numclass (length ni))
        (covw (make-array numclass))
        (covwsqrt (make-array numclass))
        (cw (make-array numclass))
        (aw (make-array numclass))
        (resw (make-array numclass))
        (vw (make-array numclass))
        (vkw (make-array numclass))
        (hw (make-array numclass))
        (iter 0) (iter1 0) (rss0 0) (rss 0) 
        (g0 0) (g1 0) (g2 0)
        (initb 0) (newb 0) b 
    )
    
;;; initialization.         
  (setf initb (bind-rows (repeat 1 p)))
  (do ((j 1 (+ j 1))) ((> j (- d 1)))
    (setf b (repeat 0 p))
    (setf (nth j b) 1)
    (setf initb (bind-rows initb b)))
  (setf rss0 0)
  (dotimes (c numclass)
    (setf v (send self :class-slice-means c))
    (setf cov (first (send self :class-cov-matrix c)))
    (setf slicesw (send self :class-y2-slices c))
    (setf nkw (mapcar 'length slicesw))
        
    (setf covsqrt (sqroot-matrix cov))    
    (setf (select covw c) cov)
    (setf (select covwsqrt c) covsqrt)
    (setf (select vw c) (matmult (diagonal (sqrt nkw)) v (inverse covsqrt)))
    (setf (select aw c) (matmult  initb covsqrt))
    (setf proj (matmult (select vw c) (transpose (select vw c))))
    (setf rss0 (+ rss0 (sum (diagonal proj))))
  )

    ;;; if d<=0, out.
    (if (< d 0) (error "d negative.") )
    (if (> d 0)
     (tagbody 
    ;;; preparation. initialize resw, g0, and iter.

  (if (= d 1)
      (tagbody
        (setf g0 0)
        (dotimes (c numclass)
            (setf v1 (select vw c)) 
            (setf a1 (select aw c))
            (setf c1 (/ (matmult v1 (transpose a1)) (sum (* a1 a1)))) 
            (setf res1 (- (transpose v1) (outer-product a1 c1)))
            (setf (select resw c) res1)
            (setf (select cw c) c1)
            (setf g0 (+ g0 (sum (^ res1 2))))
        )            
      )
      (tagbody
          (setf g0 0)
          (dotimes (c numclass)
            (setf v1 (select vw c))
            (setf a1 (select aw c))
            (setf c1 (matmult (inverse (matmult a1 (transpose a1)))
             a1 (transpose v1)))
            (setf res1 (- (transpose v1) (matmult (transpose a1) c1)))
            (setf (select resw c) res1)
            (setf (select cw c) c1)
            (setf g0 (+ g0 (sum (^ res1 2))))
        )
      )
    )
  (setf iter 0)
 
 ;;; two-layer loop.   
    (do ((i 1 (+ i 1))) ((> i itmax))
    
    ;;; begin inner loop. which updates all d rows of initb once.;;;
     
     (do ((j 1 (+ j 1))) ((> j d))
       (if (= d 1)
        (tagbody
         (dotimes (c numclass)
            (setf (select hw c) (select  cw c))
            (setf (select vkw c) (transpose (select vw c)))
          )
         )
        (tagbody
         (dotimes (c numclass)
          (setf (select hw c)(matmult (combine 1 (repeat 0 (- d 1)))
            (select  cw c) ))
          (setf (select vkw c) (+ (select resw c) 
            (outer-product (matmult (combine 1 (repeat 0 (- d 1)))
            (select aw c)) (select hw c))))
          )
         )
       )

    ;;; calculate W1 and W2. still in inner loop.

            (setf w1 0) (setf w2 0)
        (dotimes (c numclass)
            (setf w1 (+ w1 (matmult (select covwsqrt c) (select vkw c)
            (select hw c))))
            (setf w2 (+ w2 (* (sum (* (select hw c) (select hw c)))
            (select covw c)))))

    ;;; calculate the new direction \hat{b} based on  W1 and W2.
    ;;; delete the 1st row of initb. Put \hat{b} in as the last row.
    ;;; which equivalent to update the 1st row.

            (if (= d 1) (setf b (matmult (inverse w2) w1))
                (tagbody
                    (setf initb (matmult (bind-columns (repeat 0 (- d 1))
                        (identity-matrix (- d 1))) initb))
                    (setf b (matmult (transpose initb)
                    (inverse (matmult  initb (inverse w2)
                    (transpose initb))) initb (inverse w2)))
                    (setf b (matmult (inverse w2) (- (identity-matrix
                        (select (array-dimensions w2) 0)) b) w1))
                )
            )
            (setf b (/ b (sqrt (sum (* b b)))))
            (if (= d 1) (setf initb b)
                (setf initb (bind-rows initb b)))
    ;;; update aw, cw, and resw.

  (if (= d 1)
      (tagbody
          (setf g1 0)
        (dotimes (c numclass)        
            (setf a1 (matmult (transpose initb) (select covwsqrt c)))
            (setf v1 (select vw c))
            (setf c1 (/ (matmult v1 (transpose a1)) (sum (* a1 a1))))
            (setf res1 (- (transpose v1) (outer-product a1 c1)))
            (setf (select resw c) res1)
            (setf (select cw c) c1)
            (setf (select aw c) a1)
            (setf g1 (+ g1 (sum (^ res1 2))))
        )            
      )
      (tagbody
          (setf g1 0)
      (dotimes (c numclass)
        (setf a1 (matmult initb (select covwsqrt c)))
        (setf v1 (select vw c))
        (setf c1 (matmult (inverse (matmult a1 (transpose a1)))
             a1 (transpose v1)))
        (setf res1 (- (transpose v1) (matmult (transpose a1) c1)))
        (setf (select resw c) res1)
        (setf (select cw c) c1)
        (setf (select aw c) a1)
        (setf g1 (+ g1 (sum (^ res1 2))))
       )
    )
    ))
        ;;; end of inner loop. ;;;
   
    ;;; After update one the whole matrix B.                  
    ;;; calculate RSS(g1); compare with previous result (g0). 
    
        (setf iter (+ iter 1))
        (if (< (abs (- g1 g0)) eps) (return))
        (setf g0 g1)
    )
 
 ;;; end of outer loop ;;;;
 
))
(cond ((= d 0)  (list  rss0 0 0 iter))
      ((= d 1)  (list g1 initb 
              (mapcar #'transpose (coerce cw 'list)) iter) )
      ((> d 1)  (list g1 (transpose initb) cw iter))
      )
))

;;; Preparation: define  vecc. 
(defun vecc (x)
  (check-matrix x)
  (let* ((dim1 (array-dimension x 0))
     (dim2 (array-dimension x 1)))
    (make-array (list (* dim1 dim2) 1)
    :initial-contents (transpose x))))

(defun diagonal-block (x y)
(let* (temp )
  (if (eq x nil) (setf temp y)
    (tagbody 
     (let* ( (xdim1 (first (array-dimensions x)))
       (xdim2 (second (array-dimensions x)))
       (ydim1 (first (array-dimensions y)))
       (ydim2 (second (array-dimensions y)))
       )
       (setf x (bind-rows x (make-array (list ydim1 xdim2) :initial-element 0)))
       (setf y (bind-rows (make-array (list xdim1 ydim2) :initial-element 0) y))
       (setf temp (bind-columns x y))
     )))
  temp
))

(defun sqroot-matrix (x)
   (check-matrix x)
   (if (not (=  x (transpose x))) (error ("Not a symmetric matrix!")))
   (let* ((temp (eigen x)) 
      (d (nth 0  (array-dimensions x)))
      (dmat (diagonal (sqrt (first temp))))
      (vmat (make-array (list d d) 
                :initial-contents (second temp))))
   (matmult (transpose vmat) dmat vmat)
))

;;;
;;; Data set proto
;;;

(defmeth dataset-proto :fit-gpartial-sir ()
  (let* ((m1 3)
         (m2 3)
         (m3 3) 
         (name (send self :get-model-name "GP.sir"))

         (old-mod (send self :last-model :type :gpartial-inverse))
        
         (names (send self :sir-namelist))

         (n (if old-mod (reverse (set-difference names
                 (combine (send old-mod :predictors)
                          (send old-mod :response1-name)
                          (send old-mod :response2-name)))) names))
         (om (if old-mod (send old-mod :args) nil))
         
         (ans (apply #'gpartial-sir-dialog names :name name 
                :candidates n  :nslice1 m1 :nslice2 m2 :gpdim m3 om)))
    (if ans (apply #'send self :make-gpartial-sir ans))))

;;;
;;; Make GPartial SIR
;;;

(defmeth dataset-proto :make-gpartial-sir
  (&rest args &key predictors response1 response2 weights 
    nslice1 nslice2 gpdim name)
  (let* ((x (apply #'append (mapcar #'(lambda (n) (send self :data n))
                predictors)))        
         (predictor-names (combine (mapcar #'(lambda (a) (send self :labels a))
                  predictors)))
         (y1 (first (send self :data response1)))
         (response1-name response1)
         (y2 (first (send self :data response2)))
         (response2-name response2)
         (weight-name weights)
         (weights (if weights (first (send self :data weights))))
         (dataset self)     
         (r (send gpartial-sir-proto :new
                  :x (apply #'bind-columns x)
                  :y1 y1
                  :y2 y2
                  :weights weights
                  :method-response :y
                  :nslice1 nslice1
                  :nslice2 nslice2
                  :gpdim gpdim
                  :predictor-names predictor-names
                  :response1-name response1-name
                  :response2-name response2-name
                  :weights-name weight-name
                  :data self)))
    (set (intern (string-upcase name)) r)
    (send r :type :gpartial-inverse)
    (send r :args args)
    (send r :delete-slot 'included)
    (defmeth r :included (&rest args) (apply #'send dataset :included args))
    (send r :menu name)
    (when *show-display-fit* (send r :display))
    (send self :models r)
    r))
    
(defmeth gpartial-sir-proto :display ()
  (flet
    ((p-num (a b) (format t "~8,3f~8,3f " a b))
     (p-eval (a)  (format t "~11,3f" a))
     (p-espace () (format t "     "))
     (p-dir (j)  (format t "    Lin Comb ~a   " j))
     (p-head (j) (format t "    Raw    Std.  "))
     (new-line () (format t "~%"))
     (normalize (x) (coerce (/ x (sqrt (inner-product x x))) 'list)))
    (when (send self :data) (send self :description1))
    (send self :description2)
    (let* (
           (dim (send self :gpdim))
           (evectors (column-list (send self :gpdir)))
           (sds (mapcar #'standard-deviation  
                (column-list (send self :sir-x))))
           (std-dir (mapcar #'(lambda (v) 
                (normalize (* v sds))) evectors))    
           (dir (mapcar #'normalize evectors))
            (num (iseq dim))
      (names (send self :predictor-names))
      (nw (1+ (max (mapcar #'length names) 12))))
      (format t "Std. coef. use predictors scaled to have SD equal to one.~%")
      (format t (format nil "~~~dt" nw) )
      (mapcar #'p-dir (iseq 1 dim)) (new-line)
      (format t "Predictors")
      (format t (format nil "~~~dt" nw) )
      (mapcar #'p-head (iseq 1 dim)) (new-line)
      (mapcar #'(lambda (name dir sdr)
        (format t (format nil "~~~da" (1- nw)) name)
        (mapcar #'p-num (select dir (iseq dim)) (select sdr (iseq dim)))
        (new-line))
         names (transpose dir) (transpose std-dir))
      (new-line)
      (send self :ir-print-angles dir dim nw)
      (send self :test)    
    )))

;;;
;;; Description1 and Decription 2 overwrites
;;;

(defmeth gpartial-sir-proto :description1 ()
  (let ((notinc (which (mapcar #'null (send self :included)))))
    (format t "~%General Partial SIR Regression")
    (format t "~%Name of Dataset = ~a" (send (send self :data) :name))
    (format t "~%Name of Fit = ~a" (send self :name))
   (format t "~%Response (Y) = ~a" (send self :response2-name))
   (format t "~%Conditioning variable (W) = ~a" (send self :response1-name))
        
    (format t "~%Predictors = ~a" (send self :predictor-names))
    (if (send self :weights)
        (format t "~%Weights = ~a" (send self :weights-name))
        (format t "~%"))
    ;(format t "Response = ~a~%"(send self :method-response))
    (when notinc  
       (cond ((<= (length notinc) *print-n-deleted-cases*)
              (format t "Deleted cases are:~%~a~%" 
                     (select (send self :case-labels) notinc)))
             (t (format t "~a cases have been deleted.~%" (length notinc)))))))


(defmeth gpartial-sir-proto :description2 ()
  (format t "~%Slices for W     = ~a" (send self :actual-nslices1))
  (format t "~%Slices for Y within W   = ~a" (send self :actual-nslices2))
  (format t "~%Total slices     = ~a~%" (send self :actual-nslices))
  )

;;;
;;;  Plots
;;;

(defmeth gpartial-sir-proto :direction (kth &key (add nil))
"Method args: (kth)
Returns the k-th direction (linear combination) from GP.SIR."
  (let* ((f (send self :sir-x :centered t :all t))
         (evecs (column-list (send self :gpdir)))
         (name (send self :name))
         (lname (+ 1 kth))
         name dir datalist)
    (if (>= kth (length evecs)) (return-from :direction))

    (setf dir (coerce (matmult f (select evecs kth)) 'list))
    (setf name (format nil "~a.p~a" (send self :name) lname))
    (setf (slot-value 'dir-names) (cons name (slot-value 'dir-names)))
    (cond (add (setf datalist (send datalist-proto :new ':variate
                  :data dir :name name
                  :info (format nil "~a predictor ~a"
                              (send self :method) lname)))
               (send (send self :data) :add-datalist datalist))
          (t dir))))


(defmeth gpartial-sir-proto :plot-directions ()
"Method args: ()
Return the spin-plot of GP.SIR."
;;;  (if (send self :needs-computing) (send self :compute))
  (let* (
    (title (send self :title))
    (dim (send self :gpdim))
    (plotd (if (= dim 1) 1 2))
    (names (if (= plotd 1) 
        (list "Dir1" (string-downcase (send self :method-response)))
        (list "Dir1" (string-downcase (send self :method-response)) "Dir2")))
    (data #'(lambda ()
        (let* (
            (f (send self :sir-x :center t :all t))
            (evecs (column-list (send self :gpdir)))
            (dir (mapcar #'(lambda (k) (coerce (matmult f k) 'list))
                            (select evecs (iseq plotd))))
            (y (send self :sir-y :all t)))
            (append (list (first dir) y) (rest dir)))))
    (graph (send self :make-plot data names 
        :title (format nil "(~a) ~a" (send self :name) title)
        :plot-control t)))
    graph))