/* SAS program for UVA RCSC SAS/IML talk, Mar 23, 2005
 * BY: Kathy Gerber 
 * An assortment of IML commands accompanying the talk
 * to appear at 
 * http://www.itc.virginia.edu/research/talks/sasiml.html
 * 
 * ASSUMES: -- SAS, SAS/IML are available
 *          -- data set gpa.sas7bdat is in library1   
 *          -- Windows environment - minor changes
 *             should adapt program to Unix variants
 *
 * TO DO:  Set library1 to the path containing the
 *         dataset gpa
*/

%let library1 = "C:\mydat\WORK\IML\WorkShop\Data Sets";

PROC IML ;  

*--Matrix names can be from 1 to 8 characters long -- ; 
  a={1 7 3, 
      6 2 4};
   b={0 8 2, 4 1 3};
   s ={1 2 3, 4 0 6, 7 8 9} ;
  scalarB = 7.1 ;
  answer={[2]  'Yes', [2]  'No'};

*--  Operators ---;
   c0=a`;         *--- transpose    --- ;
   
   y=INV(s);     *-- Matrix inverse -- ;
   c = a*c0 ;    *-- Matrix multiplication      -- ;
   c1 = a#a ;    *-- Elementwise multiplication -- ;
   c2 = a||a ;   *-- Horizontal concatenation   -- ;
   c3 = a//a ;   *-- Vertical concatenation     -- ;
   c4 = a <> b ; *-- Choose maximimum ------------ ;
   c5 = a > b  ; 
   c6 = shape(a,6,2) ;
   c7 = shape({'Fill' 'up' 'this' 'matrix!'},2,2) ;
*----- System Action Commands --------------------------- ;
;
PRINT a, b, c0, c, c1, c2, c3, c4, c5, c6, c7, y, scalarB, answer;
FREE  c1 c2 c3 c4 c5 c6 c7 ;
SHOW a ;

*-- default storage catalog is WORK.IMLSTOR --- ;
SHOW STORAGE ;
LIBNAME mylib &library1 ;
RESET STORAGE=mylib.mystor;
STORE a ;    *--- storing a matrix -- ;
SHOW STORAGE ;
*--- Define a module ----- ;
START mymod(a);
 a = 2*a;
 PRINT a;
FINISH;

CALL mymod(a) ;

*--- Storing a module --- ;
STORE MODULE=mymod ;

*-- Nesting Modules ----- ;
t = 0 ;
START a; 
      start b; 
         t=t+1; 
      FINISH b; 
      RUN b; 
   FINISH a; 
   RUN a;
   PRINT t ;

t = 99 ;
START b; 
      t=t+1; 
   FINISH ; 
   START a; 
   RUN b; 
   FINISH; 
   RUN a;
PRINT t;

*------ A FEW FUNCTIONS    -------- ;
*------ See online documentation -- ;
*------ For a comprehensive list -- ;

* Matrix inquiry functions ;
c = NCOL(a);
PRINT c;

SHOW answer;

* Scalar functions ;
  c = log(a);
  x=SQRT(a);
  PRINT c, x ;


*-- DATA MANAGEMENT COMMANDS ----- ;

USE mylib.gpa ;
LIST all ;    *-- lists the data set --- ;
LIST POINT 12 ;
LIST NEXT ;
LIST POINT (14:22);
LIST var{gpa, hsm, hss} POINT (14:22) WHERE(GPA > 4.9);

READ ALL VAR{hsm hss} into mymat ;
PRINT MYMAT ;

*-- Create temporary dataset mydata ---- ;
varnames = {'h11', 'h22'} ; 
   CREATE mydata FROM mymat [COLNAME=varnames] ; 
   APPEND from mymat ;
CLOSE mydata ;

*-- default storage catalog is WORK.IMLSTOR --- ;

SHOW STORAGE ;
RESET STORAGE=work.mystor;
STORE MODULE=mymod ;
STORE a ;
RESET STORAGE = mylib.mystor ;


*----- apply ------------- ;

start abc(x); 
      r=x+100; 
      return (r); 
   finish abc; 

/* Modules are stored in the form of their compiled code. 
*/
 
STORE MODULE = abc ;
a={6 7 8, 
      9 10 11}; 
STORE a ;

*---- To explore persistence --- ;
*--   across SAS sessions   ---- ;
*---- RESTART SAS Session ------ ;
*---- and begin here ----------- ;

%let library1 = "C:\mydat\WORK\IML\WorkShop\Data Sets";
LIBNAME mylib &library1 ;
PROC IML ;
RESET STORAGE=mylib.mystor;
*-- Modules are in parentheses --- ;
*-- Matrices are listed next ----- ;
LOAD module=(abc) a ;



   m=apply("ABC",a);
   PRINT  m ;
   b = 100 ;
   m1=apply("abc",b);
   PRINT m1 ;

**** QUIT;

*-- CONTROL STATEMENTS ---- ;

 y = 2 ; z = 0 ;
 if y>z then 
      do; 
         if z=0 then 
            do; 
               z=4; 
               x=2#y; 
            end; 
      end;
	  
PRINT x,y,z;

 START MOD3 ;
 MSG ="Enter RESUME to continue;"; 
  PAUSE msg;
  FINISH ;
  CALL MOD3 ;
*--- The message appears in the Output Window ---- ;
*--- Highlight and run RESUME in the Editor ------ ;
  RESUME ;

*-- Running a data step or proc --- ;
*-- terminates iml automatically -- ;


DATA newdat ;
INPUT b c @@ ;
DATALINES;
2 4 3 5 6 8
5 6 7 8 3 6
;
run ;

PROC MEANS DATA = newdat  n ;
  VAR b ;
RUN ;


PROC IML ;
START SYSTEM(COMMAND); 
      CALL PUSH(" x '",command,"'; resume;"); 
      PAUSE; 
   FINISH; 
   RUN SYSTEM('cd c:\');
   RUN SYSTEM('dir');

***QUIT ;

*-- Graphics subsystem in IML ------- ;
CALL GSTART;                                                                                                                 
                                                                                                                             
x = UNIFORM(i(5)) # normal(i(5)) * 1000;       * random x data;                                                              
y = UNIFORM(i(5)) # normal(i(5)) * 1000;       * random y data;                                                              
                                                                                                                             
RUN GXYPLOT (x,y);                                                                                                           
                     
QUIT ;





****  From SAS for Monte Carlo Studies: A Guide for Quantitative Researchers   ***** ;
****  By Xitao Fan, Akos Felsovalyi, Stephen A. Sivo, and Sean C. Keenan       ***** ;
****  Program 4.9:  Generate Sample of Non-Normal Multivariate Data            ***** ;

DATA A (TYPE=CORR);  _TYPE_='CORR';  
  INPUT X1-X3;
CARDS;
1.0000  .    .
 .7402 1.0000  .
 .2054  .4173 1.0000
;

    * obtain factor pattern matrix for later data generation;

PROC FACTOR N=3 OUTSTAT=FACOUT;
DATA PATTERN; SET FACOUT;
  IF _TYPE_='PATTERN';
  DROP _TYPE_ _NAME_;
RUN;

PROC IML;                          
   USE PATTERN;                   * read in the factor pattern as a matrix 'F';
   READ ALL VAR _NUM_ INTO F;
  F=F`;           

DATA=RANNOR(J(10000,3,0));    *** generate data matrix (10000´3);
DATA=DATA`;                   *** transpose data matrix (3´10000);
Z = F*DATA;                   *** impose inter-correlations;
Z = Z`;                       *** transpose data matrix back (10000´3);

                                        * Fleishman non-normality transformation;

X1 = -.124833577 + .978350485*Z[,1] + .124833577*Z[,1]##2 + .001976943*Z[,1]##3;
X2 =  .124833577 + .978350485*Z[,2] - .124833577*Z[,2]##2 + .001976943*Z[,2]##3;
X3 = -.096435287 + .843688891*Z[,3] + .096435287*Z[,3]##2 + .046773413*Z[,3]##3;

X1=X1*15 + 100;                         * linear transformation for mean & std;
X2=X2*10 + 50;
X3=X3;
Z=X1||X2||X3;
CREATE A FROM Z [COLNAME={X1 X2 X3}];   * output a SAS working data 'A';
APPEND FROM Z;

                                        * obtaining descriptive stats for sample data;
PROC MEANS DATA=A N MEAN STD SKEWNESS KURTOSIS;
  VAR X1 X2 X3;
PROC CORR DATA=A NOSIMPLE NOPROB;
  VAR X1 X2 X3;
RUN; QUIT;
*******************************************************************;

PROC UNIVARIATE DATA = a ;
  HISTOGRAM x1 x2 x3 ;
RUN ;


*-- Using ODS with PROC IML ---------- ;
PROC DATASETS LIB =work nodetails nolist; 
DELETE myTest; 
RUN ; 
QUIT ;

ODS HTML ;
PROC IML ; 
START main; 
  myTest = {45 76, 54 23, 26 42}; 
  PRINT myTest; 
FINISH main;
 
RUN main; 
QUIT ;
 
ODS HTML CLOSE ;

********************************************** ;


/*SAS/IML algorithm for generating all n! permutations of n objects (Tobias, 1995). */

proc iml;
start perm(n);
	avail = shape(1,1,n);
 	idx   = shape(0,1,n);
	np = 1;
	do i = 2 to n; np = np * i; end;
		perm = shape(0,np,n);
      		ip = 0;
      		lev = 1;
      		do while (lev);
         		sym = loc(avail & ((1:ncol(avail)) > idx[lev]));
		if (ncol(sym)) then do;
            		idx[lev] = sym[1];
            		avail[idx[lev]] = avail[idx[lev]] - 1;
            		if (lev < n) then do;
               			lev = lev + 1;
            		end;
			else do;
               			ip = ip + 1;
               			perm[ip,] = idx;
               			lev = lev - 1;
               			avail[idx[lev:(lev+1)]] = avail[idx[lev:(lev+1)]] + 1;
               			idx[lev+1] = 0;
			end;
		end;
         		else do;
			lev = lev - 1;
            		if (lev) then do;
               			avail[idx[lev]] = avail[idx[lev]] + 1;
               			idx[lev+1] = 0;
               		end;
            	end;
	end;

	return(perm);
finish;

x=perm(4);	/* Generates all 4! = 24 permutations of the set {1,2,3,4} */
print x;
quit;