/* Program to accompany the article
   Wicklin, R. (2016) "Visualize the Cantor function" published 29JUN2016
   on The DO Loop blog: 
   http://blogs.sas.com/content/iml/2016/06/29/viz-cantor-function-in-sas.htm
*/
proc iml;
/* rows of B contain all 8-digit combinations of 0s and 2s */ 
B = expandgrid({0 2}, {0 2}, {0 2}, {0 2}, {0 2}, {0 2}, {0 2}, {0 2});
B = B[2:nrow(B),];      /* remove first row of zeros */
k = ncol(B);            /* k = 8 */
 
v = 3##(-(1:k));       /* vector of powers 3^{-i} */
t = B * v`;            /* x values: right-hand endpts of middle-third intervals */
 
u = 2##(-(1:k));       /* vector of powers 2^{-i} */
f = B/2 * u`;          /* y values: Cantor function on Cantor set */
call series(t, f);     /* graph the Cantor function */
quit;


/* More elaborate program that adds tick marks to the graph */
proc iml;
B = expandgrid({0 2}, {0 2}, {0 2}, {0 2}, 
               {0 2}, {0 2}, {0 2}, {0 2});
*B = expandgrid({0 2}, {0 2}, {0 2}, {0 2});
B = B[2:nrow(B),];  /* remove first row of zeros */
k = ncol(B); 
v = 3##(-(1:k));  /* vector of powers 3^{-i} */
t = B * v`;       /* some points in the Cantor set */

u = 2##(-(1:k));  /* vector of powers 2^{-i} */
f = B/2 * u`;     /* values of Cantor function on Cantor set */
title "Cantor Function";
call series(t, f) yvalues=do(0,1.01,0.125);

/* For each row of a matrix A, find the rightmost column that contains
   a positive number. For example:
   A = {1 0 1 1 0 0 0,     
        0 1 0 1 1 1 0,     
        1 0 0 0 0 0 0};    
   p = RightmostPosCol(A); * = {4, 6, 1} ;
   The reason:
   - In the first row, column 4 is the last nonpositive value in row
   - In the second row, column 6 is the last nonpositive value in row
   - In the third row, column 1 is the last nonpositive value in row
*/
start RightmostPosCol(A);
   cols = j(nrow(A), 1, 0);
   do j = 1 to ncol(A);
      idx = loc(A[,j] > 0);
      if ncol(idx)>0 then cols[idx] = j;
   end;
   return (cols);
finish;

/* remove first row of zeros */
depth  = RightmostPosCol(B);
z = t - 3##(-depth); /* left hand endpts of middle-third intervals */
x = 0 // colvec(z || t) // 1;    /* points in Cantor set */
y = 0 // colvec(f || f) // 1;    /* function values */

/* extract labels for the first few intervals */
d = 1 // colvec(depth || depth) // 1; /* step at which interval was added */
idx = loc(d <= 2);
xvalues = x[idx];
xlabels = '"' + putn(xvalues, "fract5.") + '"';
idx = loc(d <= 3);
yvalues = unique(y[idx]);
ylabels = '"' + putn(yvalues, "fract5.") + '"';

*print xvalues xlabels yvalues ylabels;
create Cantor var {x y}; append; close;

submit xvalues xlabels yvalues ylabels;
proc sgplot data=Cantor2;
series x=x y=y;
xaxis values=(&xvalues) valuesdisplay=(&xlabels) grid;
yaxis values=(&yvalues) valuesdisplay=(&ylabels) grid;
run;
endsubmit;
quit;