/* Program to accompany the blog post:
   "Mathematical art: Weaving matrices,"
   by Rick Wicklin, published Sep 9, 2015,
   The DO Loop blog, 
   http://blogs.sas.com/content/iml/2015/09/02/weaving.html
*/

/* Heat maps use gray for missing values. Define
   style that uses the background color (white) */
proc template;
   define style mystyle;
      parent = styles.htmlblue;
      class GraphMissing / color = GraphColors('gwalls');
   end;
run;
ods html style=mystyle;
ods graphics / width=300px height=300px;

proc iml;
/* Experiment with the weaving algorithm on numeric matrix */
reset linesize=128;
n = 10;
v = "1":char(n);
x = shape(1:n##2, n); /* original matrix */

w = j(n, n, .);    /* start with all missing values */
/* first half of columns ==> vertical strips */
j = 1:n/2;         
w[, 2*j] = x[, j]; /* copy to every other column */
print w[c=v r=v];

/* second half of columns ==> horizontal strips */
i = (n/2+1):n;
w[2*j,] = x[, i]`; /* copy to every other row */
print w[c=v r=v];

/* weaving: raise certain rows from vertical strips */
k = do(1, n/2, 2); 
w[2*k, 2*k] = x[2*k, k]; /* every other for even rows */
k = do(2, n/2, 2);
w[2*k, 2*k] = x[2*k, k]; /* every other for odd rows */
print w[c=v r=v];
quit;

proc iml;
/* encapsulate the weaving algorithm into a module */
start weave(x);
   n = nrow(x);         /* n must be even */
   y = j(n, n, .);
   j = 1:n/2;           /* vertical strips */
   y[, 2*j] = x[, j];   /* copy to every other column */
   i = (n/2+1):n;       /* horizontal strips */
   y[2*i-n,] = x[, i]`; /* copy to every other row */
   k = do(1, n/2, 2);   /* weaving: raise certain rows */
   y[2*k, 2*k] = x[2*k, k]; /* every other for even rows */
   k = do(2, n/2, 2);
   y[2*k, 2*k] = x[2*k,k];  /* every other for odd rows */
   return( y );
finish;

x = shape(1:100, 10); /* original 10x10 matrix */
w = weave(x);
print w;


/* helper module: rotate matrix 90 degrees counterclockwise */
/* http://blogs.sas.com/content/iml/2013/10/18/rotating-matrices.html */
start Rot90(m);
   return( T(m[,ncol(m):1]) );         /* left-right flip, then transpose */
finish;

/* now assign values (colors) to original matrix and see what happens */
/* some alternative color ramps are below: */
/*
ramp = "rainbow";
ramp = "temperature";
ramp = {CX1B485B CX342D1C CX532414 CX764B32 CXDC7B10 CXD4AB56 CX465D4C CX83A089};
*/
/* use muted "rainbow" ramp: {magenta cyan yellow red} */
ramp = {CXAA00AA CX00AAAA CXAAAA00 CXAA0000};

/* some options for every heat map */
%let opt = colorramp=ramp displayoutlines=0 showlegend=0;

/* create 20x20 grid of points in [-1,1]x[-1,1] */
n = 20;
t = do(-1, 1, 2/(n-1));
xy = expandgrid( t, t );
x = xy[,1];  y = xy[,2];

/* Create mathematical surfaces ("contours") and weave */
/* linear function: z = x */
z = x;
m = Rot90(shape(z, n));
call heatmapcont(m) title="Original painting = 'lines'" &opt;
w = weave(m);
call heatmapcont(w) title="Weaving of 'lines'" &opt;


/* linear function: z = x + y */
z = x + y;
m = Rot90(shape(z, n));
call heatmapcont(m) title="Original painting = 'linear'" &opt;
w = weave(m);
call heatmapcont(w) title="Weaving of 'linear'" &opt;


/* quadratic bowl: z = x##2 + y##2 */
z = x##2 + y##2;
m = Rot90(shape(z, n));
call heatmapcont(m) title="Original painting = 'bowl'" &opt;
w = weave(m);
call heatmapcont(w) title="Weaving of 'bowl'" &opt;


/* quadratic saddle: z = x##2 - y##2 */
/* This is the SECRET weave! */
z = x##2 - y##2;
m = Rot90(shape(z, n));
call heatmapcont(m) title="Original painting = 'saddle'" &opt;
w = weave(m);
call heatmapcont(w) title="Weaving of 'saddle'" &opt;


/* Concentric rings: z = sin( pi*(x##2 + y##2) ) */
pi = constant('pi');
z = sin(pi*(xy[,1]##2 + xy[,2]##2));
m = Rot90(shape(z, n));
call heatmapcont(m) title="Original painting = 'rings'" &opt;
w = weave(m);
call heatmapcont(w) title="Weaving of 'rings'" &opt;

/****************************************************************/
/* You can play with the color ramps, the number of strips, etc */
/* try different color ramp */
ramp = {CX1B485B CX342D1C CX532414 CX764B32 CXDC7B10 CXD4AB56 CX465D4C CX83A089};
z = x##2 + y##2;
m = Rot90(shape(z, n));
w = weave(m);
call heatmapcont(w) &opt;

quit;