/* Program to accompany the blog post: "Mathematical art (part 2): Unweaving matrices," by Rick Wicklin, published Sep 4, 2015, The DO Loop blog, http://blogs.sas.com/content/iml/2015/09/04/unweave.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; /* 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; /* Because of the weaving, we have lost information. Some values are lost forever. */ start unweave(w); n = nrow(w); y = j(n,n,.); j = 1:n/2; y[,j] = w[, 2*j]; /* compress the vertical strips */ i = (n/2+1):n; y[,i] = w[2*j,]`; /* rotate and compress the horizontal strips */ /* Because of the weaving, some values are unknown. Insert missing. */ k4 = do(2, n, 4); k2 = do(2, n/2, 2); y[k4, k2] = .; /* 1st half of even rows that are not divisible by 4 */ k2 = do(n/2+1, n, 2); y[k4, k2] = .; /* 2nd half of even rows that are not divisible by 4 */ k4 = do(4, n, 4); k2 = do(1, n/2, 2); y[k4, k2] = .; /* 1st half of even rows divisible by 4 */ k2 = do(n/2+2, n, 2); y[k4, k2] = .; /* 2nd half of even rows divisible by 4 */ return ( y ); finish; /* test the unweaving algorithm on a 10x10 ordered matrix */ x = shape(1:100, 10); /* original 10x10 matrix */ w = weave(x); print w; z = unweave(w); print z; /* 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 */ /* 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]; /* 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 range={-1 1}; w = weave(m); call heatmapcont(w) title="Secret Weave" &opt range={-1 1}; /* unweave the secret matrix */ u = unweave(w); call heatmapcont(u) title="Unweaving" &opt range={-1 1}; /* 2-D interpolation by piecewise constant function. Assume f(x,y) is continuous and that A[i,j] = f(x_i, y_j). Find missing values in A and interpolate based on average values of nearest neighbors. */ A = u; /* specify location of missing values */ LL = n*(n-1)+1; /* lower left corner*/ LR = n*n; /* lower right corner*/ R = do(4*n, n##2-1, 4*n); /* right edge */ L = R - n + 1; /* left edge */ B = do(LL+2, LL + n/2 -1, 2) || do(LL + n/2 +1, n*n-2, 2); /* bottom edge */ /* 1. interpolate corners */ nbrs = (LL-n) // (LL+1); A[LL] = mean( A[nbrs] ); nbrs = (LR-n) // (LR-1); A[LR] = mean( A[nbrs] ); /* 2. interpolate left edge */ do i = 1 to ncol(L); nbrs = (L[i]-n) // (L[i]+1) // (L[i]+n); A[L[i]] = mean( A[nbrs] ); end; /* 3. interpolate right edge */ do i = 1 to ncol(R); nbrs = (R[i]-n) // (R[i]-1) // (R[i]+n); A[R[i]] = mean( A[nbrs] ); end; /* 4. interpolate bottom edge */ do i = 1 to ncol(B); nbrs = (B[i]-1) // (B[i]-n) // (B[i]+1); A[B[i]] = mean( A[nbrs] ); end; /* 5. interpolate interior */ Intr = loc(A=.); do i = 1 to ncol(Intr); nbrs = (Intr[i]-1) || (Intr[i]+1) || (Intr[i]-n) || (Intr[i]+n); A[Intr[i]] = mean( A[nbrs] ); end; /* show the final estimate for the original artwork */ call heatmapcont(A) title="Interpolation" &opt range={-1 1}; /* for comparison, show the original artwork */ call heatmapcont(m) title="Original image" &opt range={-1 1}; /* by how much does the approximation differ from the truth? */ Diff = A - m; maxDiff = max(abs(Diff)); print maxDiff; RelDiff = Diff / m; maxRelDiff = max(abs(RelDiff)); print maxRelDiff; MeanSqErr = mean( colvec(Diff##2) ); print MeanSqErr; /* plot location of differences, mostly on edges */ call heatmapcont(RelDiff) title="Difference matrix" displayoutlines=0 ;