/* 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;