/* Program to accompany the article "The dunkard's walk in 2-D,"
   Rick Wicklin, published 12AUG2015, The DO Loop blog,
   http://blogs.sas.com/content/iml/2015/08/05/2d-drunkards-walk.html
*/
proc iml;
call randseed(54321);
dim = 2;                   /* number of directions */
NumSteps = 10;             /* number of steps */

coord = j(NumSteps,1,.);   /* allocate vector for random direction */
call randgen(coord, "Table", j(1,dim,1)/dim); /* random directions */
incr = sample({-1 1}, 1 // NumSteps);         /* random step +/-1 */
step = T(1:NumSteps);
print step coord incr;

x = j(NumSteps, dim, 0);      /* (N x d) array */
idx = sub2ndx(dimension(x), step||coord);
x[ idx ] = incr;
print x;

/***************************************************/
proc iml;
start vcusum(x);
   c = j(nrow(x), ncol(x), .);
   do i = 1 to ncol(x);
      c[,i] = cusum(x[,i]);
   end;
   return (c);
finish;

start DrunkardWalk(dim, NumSteps);
   coord = j(NumSteps,1,.);   /* allocate vector for random direction */
   call randgen(coord, "Table", j(1,dim,1)/dim); /* random directions */
   incr = sample({-1 1}, 1 // NumSteps);         /* random step +/-1 */
   step = T(1:NumSteps);
   x = j(NumSteps, dim, 0);      /* (N x d) array */
   idx = sub2ndx(dimension(x), step||coord);
   x[ idx ] = incr;
   return( x );
finish;
store module=_all_;

call randseed(54321);
dim = 2;                   /* number of directions */
NumSteps = 10;             /* number of steps */
x = DrunkardWalk(dim, NumSteps);
Y = vcusum(j(1,dim,0) // x);                  /* show initial position at t=0 */
Q = t(0:nrow(x)) || Y;
create walk from Q[c={"Step" "x" "y"}];
append from Q;
close;
quit;


title "A Single Drunkard's Walk";
proc sgplot data=walk;
series x=x y=y / datalabel=Step markers;
refline 0/ axis=x;
refline 0/ axis=y;
xaxis integer min=-2 max=2; 
yaxis integer min=-2 max=2; 
run;




/************************************************/
proc iml;
load module=(vcusum DrunkardWalk);
dim = 2;
NumSteps = 6;

call randseed(1234);
NumWalks = 10000;
finalPos = j(NumWalks, dim);
do i = 1 to NumWalks;
   x = DrunkardWalk(dim, NumSteps);
   finalPos[i,] = x[+,];
end;

/* summarize and visualize 2-D data */
create walk from finalPos[c={"x" "y"}];
append from finalPos;
close;

submit;
   proc freq data=walk noprint;   /* compute frequency distribution */
   tables x*y / list sparse out=freqOut;
   run;
endsubmit;

/* read into IML vectors */
use freqOut;
read all var {x y count};
close;

nX = range(x) + 1;
nY = range(y) + 1;
W = shape(count, nX, nY);   /* reshape into rectangular matrix */

ods graphics / width=600px height=500px;
title "10,000 Drunkard's Walks";
title2 "Distribution of Counts for Final Positions";
call HeatmapCont(W) xValues = unique(x) yValues=unique(y);

pct = W / sum(W);
print pct;
quit;