;
; singular value decomposition exercise.  This IDL program illustrates
; the math involved creating the inverse Interaction matrix.
;
; For AO, we can determine an Interaction matrix (Im), so that:
;
;   Sv[s] = Im[s][a] * Av[a]
;
; The Im is determined by execising each actuator and measuring
; the sensor response. The would like to invert Im to gave a
; control matrix Iim, so that,
;
;   Av[a] = Iim[a][s] * Sv[s]
;
; To determine Iim, we use singular value decomposition.
; svd(Im) give Um Wm Vm so that,
;
;   Im[s][a] = Um[s][a] ## Wm[a][a] ## transpose(Vm[a][a])
;
; Then we can calcuate Iim using:
;
; Iim[a][s] = Vm[a][a] ## Wim[a][a] ## Utm[a][s], where
;
;    Wim = reciprocate each diagonal element of (W)
;    Utm = transpose(Um)
;
;
PRO  svdex

   ; ------------------------------------------------
   ; make some source data
   ;
   s = 6   ; 6 sensors
   a = 6   ; 6 actuators

   ; define array I
   Im = [ [0.9544, 0.0155, 0.0027, 0.0074, 0.0014, 0.0186 ], $
          [0.0018, 0.9341, 0.0175, 0.0155, 0.0173, 0.0138 ], $
          [0.0011, 0.0087, 0.9669, 0.0099, 0.0123, 0.0011 ], $
          [0.0132, 0.0018, 0.0192, 0.9419, 0.0039, 0.0200 ], $
          [0.0091, 0.0053, 0.0117, 0.0111, 0.9548, 0.0080 ], $
          [0.0125, 0.0193, 0.0171, 0.0075, 0.0055, 0.9381 ]]

   ; ------------------------------------------------
   ; singular value decomposition:
   ;   svd(Im) give Wv, Um, Vm.
   ;   Where I = Um ## Wm ## transpose(Vm)
   ;

   print, 'Im:'
   print, Im

   SVDC, Im, Wv, Um, Vm

   print, 'Wv:'
   print, Wv

   print, 'Um:'
   print, Um

   print, 'Vm:'
   print, Vm

   ;
   ; verify, I = U * W * t(V):
   ;
   Wm = FLTARR( s, s)
   for k=0, s-1 do Wm[k,k] = Wv[k]
   R = Um ## Wm ## TRANSPOSE(Vm)
   print, 'U*W*transpose(V) (should equal I):'
   print, R

   ; ------------------------------------------------
   ;   Determine inverse I. The inverse(I) is
   ;      inverse(I) = V * inverse(W) * transpose(U)
   ;
   Wim = fltarr(s,s)
   for k=0, s-1 do begin
      if( Wm[k,k] LT 1.00000e-05 ) then begin
         Wim[k,k] = 0
      endif else begin
         Wim[k,k] = 1/Wm[k,k]
      endelse
   endfor
   Utm = transpose(Um);

   print, 'inverse(Wm)'
   print, Wim
   print, 'transpose(Um)'
   print, Utm

   Iim = Vm ## Wim ## Utm
   print, 'Iim:'
   print, Iim

   ; ------------------------------------------------
   ;  Does it work?
   ;  (Sv ## I) ## inverse(I) == Sv
   ;
   Sv = [0.48, 0.51, 0.85, 0.52, 0.56, 0.42]
   print, 'Sv:'
   print, Sv

   Av = Im ## Sv
   R = Iim ## Av
   print, 'inverse(I) * Sv'
   print, R

END