/* title 'Some useful IML modules for matrix operations'; */

 /*
    Accessing the modules:
       %include iml(matlib);

   All of these are functions, returning a value, so should be assigned
   to some matrix name, e.g.,
      u = proj(y, Z);
      r = r(A);
  
  *----- Matrix utilities -------
    row(X)    - Convert a matrix into a row vector
    col(X)    - Convert a matrix into a col vector
  *----- Det, Rank, Projections -----
    r(A)      - Returns the rank of the matrix A.
    proj(y,X) - Returns the projection of vector y on matrix X.
    minor(A,i,j) - Returns the (i,j) minor of matrix A.
    cofactor(A,i,j) - Returns the (i,j) cofactor of matrix A.
  *----- Statistical functions ------
    len(X)    - Returns the lengths of the columns of matrix X
    dev(x)    - Returns matrix X in deviations from its column means.
    scp(x)    - Returns sum of squares and cross-products(X`X).
    cov (x)   - Compute covariance matrix of columns of X  [built-in in 9.3]
    corr (X)  - Compute correlations among columns of X  [built-in in 9.3]
    median(D) - Median of each column of a matrix D
  *----- Elementary row operations ----
    rowadd (X, from, to, mult) - Add 'mult'*row 'from' to row 'to'
    rowswap(X, from, to) - Interchange rows 'from' and 'to' of matrix X
    rowmult(X, row, mult) - Multiply one 'row' by a constant 'mult'
  */

 /*-------------------------------------------------------------*
  * row(X) - Convert a matrix into a row vector
  *-------------------------------------------------------------*/
start row (x);
   if (nrow(x) = 1) then return (x);
   if (ncol(x) = 1) then return (t(x));
   n = nrow(x) * ncol(x);
   return (shape(x,1,n));
   finish;

 /*-------------------------------------------------------------*
  * col(X) - Convert a matrix into a col vector
  *-------------------------------------------------------------*/
start col (x);
   if (ncol(x) = 1) then return (x);
   if (nrow(x) = 1) then return (t(x));
   n = nrow(x) * ncol(x);
   return (shape(x,n,1));
   finish;

 /*-------------------------------------------------------------*
  * r(A) - Returns the rank of the matrix A.
  *-------------------------------------------------------------*/
start r(A);

   reset noprint;
   *-- rank = number of nonzero rows/cols in echelon form;
   e = echelon(A);
	r = (e ^= 0)[,+];    *-- rows: # of non-zero elements;
	r = (r ^= 0)[+,];    *-- # non-zero rows;
	reset print;
	return (r);

/* Two short forms:
	return ( (((echelon(A)^=0)[,+]) ^=0)[+,] );
   return ( round(trace(ginv(a)*(a))) );
*/
   finish;
 
 /*-------------------------------------------------------------*
  * minor(A,i,j) - Returns the (i,j) minor of matrix A.
  * INPUT   A   (m,m)          Square input matrix
  *-------------------------------------------------------------*/
start minor(A,i,j);
   return( det( A[remove((1:nrow(A)),i), remove((1:nrow(A)),j)] ) );
finish;
 
 
 /*-------------------------------------------------------------*
  * cofactor(A,i,j) - Returns the (i,j) cofactor of matrix A.
  * INPUT   A   (m,m)          Square input matrix
  *-------------------------------------------------------------*/
start cofactor(A,i,j);
   minor = minor(A, i,j);
   return ( ((-1)##(i+j)) # minor );
finish;
 
 /*----------------------------------------------------------------*
  * proj(y,X) - Returns the projection of vector y on matrix X.
  * INPUT   y    (n,1)         input vector
  *         X    (n,m)         input matrix
  * The function returns an (n,1) vector, the projection of y on X 
  *----------------------------------------------------------------*/
start proj(y,X) global(P);
   reset noprint;
   if nrow(X)=1 then X = t(X);
   XPX = t(X) * X;
   P = X * ginv(XPX) * t(X);
   reset print;
   return( P * y );
   finish;
 
   *-- Define a length function (LEN);
 /*-------------------------------------------------------------*
  * len(X) - Returns the lengths of the columns of matrix X
  *-------------------------------------------------------------*/
start len(X);
   return (sqrt( X[##,] ));
   finish;

 /*--------------------------------------------------------------*
  * dev(x) - Returns matrix X in deviations from its column means.
  * INPUT:   x    (n,m)         input matrix
  * OUTPUT:  d    (n,m)         output matrix
  *--------------------------------------------------------------*/
start dev(x);
*   d = x - repeat(x[:,], nrow(x), 1);
   return(x - repeat(x[:,], nrow(x), 1));
   finish dev;

 /*--------------------------------------------------------------*
  * scp(x) - Returns sum of squares and cross-products(X`X).
  * INPUT:   x    (n,m)         input matrix
  * OUTPUT:  d    (n,m)         output matrix
  *--------------------------------------------------------------*/
start scp(x);
   return(t(X) * X);
   finish;

%macro corrcov;
*-- corr() and cov() defined in SAS 9.3;
%if &sysver < 9.3 %then %do; 
 /*---------------------------------------------------------------*
  * corr (X) - Compute correlations among columns of X
  * INPUT :  X   (n x m) matrix, where n = number of data points
  *                  m = number of variables
  * OUTPUT:  returns correlation coefficients matrix  (m x m)
  *---------------------------------------------------------------*/

start corr (x);
  reset noprint;
  d = dev(x);                      * correct for means;
  xpx=t(d)*d;                      * crossproduct;
  v = vecdiag(xpx);                * diagonal values;
  v = 1/sqrt(choose(v=0,.,v));     * account for constants;
  v = choose(v=.,0,v);
  corr= diag(v) * xpx * diag(v);   * correlation matrix;
  reset print;
  return (corr);
finish;
 
 /*---------------------------------------------------------------*
  * cov (x) - Compute covariance matrix of columns of X
  *
  * INPUT :  X   (n x m) matrix, where n = number of data points
  *                  m = number of variables
  * OUTPUT:  returns covariance matrix  (m x m)
  *---------------------------------------------------------------*/
 
start cov (x);
  reset noprint;
  d = dev(x);                      * correct for means;
  xpx=t(d)*d;                        * crossproduct;
  cov = xpx / (nrow(d)-1);
  reset print;
  return (cov);
finish;
%end;
%mend;
%corrcov;
 
*-- Elementary row operations;
 /*---------------------------------------------------------------*
  * rowadd (X, from, to, mult) - Add 'mult'#row 'from' to row 'to'
  *---------------------------------------------------------------*/
start rowadd (X, from, to, mult);
	reset noprint;
   y = x;
   y[to,] = y[to,] + mult # y[from,];
   reset print;
   return(y);
   finish;

 /*---------------------------------------------------------------*
  * rowswap(X, from, to) - Interchange two rows of a matrix
  *---------------------------------------------------------------*/
start rowswap(X, from, to);
   reset noprint;
   y = x;
   temp = y[to,];
   y[to,] = y[from,];
   y[from,] = temp;
   reset print;
   return(y);
   finish;

 /*---------------------------------------------------------------*
  * rowmult(X, row, mult) - Multiply one row by a constant
  *---------------------------------------------------------------*/
start rowmult(X, row, mult);
   reset noprint;
   y = x;
   y[row,] = mult # y[row,];
   reset print;
   return(y);
   finish;

 /*---------------------------------------------------------------*
  * mpower(A, p) - Matrix power, A^p
  *---------------------------------------------------------------*/
start mpower(A, p);
   n = nrow(A);
   if nrow(A) ^= ncol(A) then do;
      end;

   if p=floor(p) then B = I(n);
   else if mod(p,1) = 0.5 then do; 
      B = half(A);
      p = floor(p);
      end;

   do i=1 to p;
      B = B * A;
      end;
   return (B);
   finish;

 /*---------------------------------------------------------------*
  * median(D) - Median of each column of a matrix D
  *---------------------------------------------------------------*/
start median(D);
   n = nrow(D);
   if n=1 then return(d);
   c = ncol(D);
   med = J(1,c,0);
   m = floor( ((n+1)/2) || ((n+2)/2) );
   do i = 1 to c;
       t = d[,i];
      rt = rank(t);
      t[rt,]=d[,i];
      med[,i] = (t[m,])[+,]/2;
      end;
   return(med);
finish;