/* Richard A. DeVenezia
* www.devenezia.com
*/
/* This code evolved from some originally posted
* to SAS-L on August 5, 2003
*
* An adjacency problem:
* ----
* Consider a member, mx, with a value in the range of A to Z.
* Consider data describing S1,...,Si,...,Sn containing
* members m{i}{1},...,m{i}{n{i}}
*
* e.g
* set id
*  1   A
*  1   B
*  1   C
*  2   D
*  3   A
*  3   F
*  4   H
*  4   I
*  4   J
*  4   K
*  5   F
*  5   E
*
* becomes
*
* superset id
* 1_3_5     A
* 1_3_5     B
* 1_3_5     C
* 2         D
* 1_3_5     A
* 1_3_5     F
* 4         H
* 4         I
* 4         J
* 4         K
* 1_3_5     F
* 1_3_5     E
*
* The problem is to find the fewest mutually exclusive
* supersets S* such that each S* contains all S sharing at
* least one member.
*
* The problem could also be construed as finding valid
* paths through linked groups.
*
* The solution demonstrates how to program a recursive
* algorithm in a SAS Data Step using GOTO and ARRAY to
* emulate functional calls with stacks.
* GOTO is required because Data Step only allows for
* 10 levels of LINK
*/

/* Step 1.
* Construct sample data containing groups in several mutually
* exclusive divisons.  Note: Four supersets is a minimum imposed
* by the construction algorithm, however there may be several
*/
data sets (keep=set id);
  length set 8 id $1 group $26;
  foo = compress
      ( '1234 1243 1324 1342 1423 1432'
     || '2134 2143 2314 2341 2413 2431'
     || '3124 3142 3214 3241 3412 3421'
     || '4123 4132 4213 4231 4312 4321'
     || '12' );
  retain seed 512;
  retain nSets 12;
  retain maxMembers 8;
  do set = 1 to nSets;
    superset = 1 + int(4*ranuni(seed));
    nMembers = 1 + maxMembers * ranuni (seed);
    group = '';
    ssc = put(superset,1.);
    do i = 1 to nMembers;
      do until (substr(foo,ix,1) = ssc);
        ix = 1 + int(26*ranuni(seed));
      end;
      substr(group,ix,1) = byte (ix+64);
    end;
    group = compress(group);
    do i = 1 to length (group);
      id = substr (group,i);
      output;
    end;
  end;
run;
/* Step 2.
* Determine number of sets in data
*/
proc sql;
  reset noprint;
  select count (distinct set) into :nSets from sets;
  %let nSets = &nSets;
quit;
/* Step 3.
* Process the data looking for supersets.
* While doing so, output information that is used to plot
* the execution state.  This will help you visualize
* what is going on.
*/
data
  supersets (keep=sset id)
  stages (keep=stage x y state)
  anno (keep=stage x1 y1 x2 y2 ix info)
;*  / debug
;
  length stage x1 y1 x2 y2 ix 8.;
  array S [ &nSets, 26 ] _temporary_;
  array row_stack [ %eval (26*&nSets) ] _temporary_;
  array col_stack [ %eval (26*&nSets) ] _temporary_;
  array claimedRow [ &nSets ] _temporary_;
  array claimedCol [ 26 ] _temporary_;
  retain ix jxix kxix 0;
  array jx_stack [ &nSets ] _temporary_;
  array kx_stack [ 26 ] _temporary_;
  /* Read all the sets data into an array
   * Mark each member as unclaimed (-1)
   */
  do while (not eog);
    set sets end=eog;
    row = set;
    col = rank (id) - 64;  * A-Z -> 1-26;
    S [ row, col ] = -1;
  end;
link puts;
  link outputStage;
  x2 = .;
  /* Examine every member in its context
   * If the member is unclaimed by a superset
   * then increment the superset and claim the
   * the row the member is in
   * Note: Claiming a row initiates a recursive
   * claim of column/row/column/...
   */
  ss = 0;
  do row = 1 to &nSets;
    if claimedRow [ row ] then continue;
    do col = 1 to 26;
      if S [ row, col ] = -1 then do;
        ss + 1;
        link claimRow;
        col = 26;
      end;
    end;
  end;
link putS;
  /* Evaluate the matrix containing the
   * superset assignments
   */
  length sset $%eval(&nSets*3);
  do ix = 1 to ss;
    sset = '';
    do i = 1 to &nSets;
      if claimedRow = ix then do;
        if sset = ''
          then sset = trim(left(put(i,4.)));
          else sset = trim(sset) || '_' || trim(left(put(i,4.)));
      end;
    end;
    do i = 1 to &nSets;
      if claimedRow = ix then do;
        do j = 1 to 26;
          if S [ i, j ] = ix then do;
            id = byte (j + 64);
            output supersets;
          end;
        end;
      end;
    end;
  end;
  stop;

outputStage:
  stage+1;
  do _row = 1 to &nSets;
  do _col = 1 to 26;
    state = S [ _row, _col ];
    y = _row;
    x = _col;
    if state then output stages;
  end;
  end;
return;
outputAnno:
  x1 = x2;
  y1 = y2;
  x2 = col;
  y2 = row;
  output anno;
return;

pushRowCol:
  ix + 1;
  row_stack [ ix ] = row;
  col_stack [ ix ] = col;
return;
popRowCol:
  row = row_stack [ ix ] ;
  col = col_stack [ ix ] ;
  ix + (-1);
return;
pushJx: jxix + 1; jx_stack [ jxix ] = jx; return;
popJx : jx = jx_stack [ jxix ]; jxix + (-1); return;
pushKx: kxix + 1; kx_stack [ kxix ] = kx; return;
popKx : kx = kx_stack [ kxix ]; kxix + (-1); return;

%*--------------------------------------------------------------------;
claimRow:
  link OutputStage; info = 'R'; x2=col; y2=row; x1=.; y1=.; output anno;
  if claimedRow [ row ] then
    goto claimRowReturn;
  claimedRow [ row ] = ss;
  col = 1;
  jx = ix;
  do until (col > 26);
    if S [ row, col ] = -1 then do;
      link pushRowCol;
      S [ row, col ] = ss;
      link outputStage;
      link outputAnno;
    end;
    col + 1;
  end;
  do while (jx < ix);
    link popRowCol;
    link pushRowCol;
    link pushJx;
    goto claimColumn;
claimColumnReturn:
    link popJx;
    link popRowCol;
  end;
  if ix = 0 then return;
  goto claimRowReturn;

%*--------------------------------------------------------------------;
claimColumn:
  link OutputStage; info = 'C'; x2=col; y2=row; x1=.; y1=.; output anno;
  if claimedCol [ col ] then
    goto claimColumnReturn;
  claimedCol [ col ] = ss;
  row = 1;
  kx = ix;
  do until (row > &nSets);
    if S [ row, col ] = -1 then do;
      link pushRowCol;
      S [ row, col ] = ss;
      link outputStage;
      link outputAnno;
    end;
    row + 1;
  end;
  do while (kx < ix);
    link popRowCol;
    link pushRowCol;
    link pushKx;
    goto claimRow;
claimRowReturn:
    link popKx;
    link popRowCol;
  end;
  goto claimColumnReturn;
%*--------------------------------------------------------------------;
putS:
  do ii = 1 to &nSets;
    put ii 2. '. ' @;
    do jj = 1 to 26;
      put S[ii,jj] 2. +1 @;
    end;
    put;
  end;
  put;
return;
format ix stage x1 y1 x2 y2 4.;
run;

/* Step 4.
* Generate gifanim that demonstrates the
* program flow performed
*/
options nosource;
*/*;
options mprint;
data annotate;
  set anno;
  retain xsys ysys '2';
  length function color $8;
  retain xx yy xx1 xx2 yy1 yy2;
  if info = 'R' and y1=. and x1=. then do;
    color='GRAY';
    xx1 =  0; yy1 = y2;
    xx2 = 27; yy2 = y2;
    xx = x2;
    yy = y2;
  end;
  else
  if info = 'C' and x1=. and y1=. then do;
    color='GRAY';
    xx1 = x2; yy1 = 0;
    xx2 = x2; yy2 = &nSets. + 1;
    xx = x2;
    yy = y2;
  end;
  x=xx-.25;y=yy-.25; function='MOVE'; output;
  x=xx+.25;y=yy+.25; function='DRAW'; output;
  x=xx-.25;y=yy+.25; function='MOVE'; output;
  x=xx+.25;y=yy-.25; function='DRAW'; output;
  x=xx1; y=yy1; function = 'MOVE'; output;
  x=xx2; y=yy2; function = 'DRAW'; size=1; output;

  if x1 & x2 & y1 & y2;
  color='BLACK';
  x = x1; y = y1; function = 'MOVE'; output;
  x = x2; y = y2; function = 'DRAW'; size=2; output;
/*
  function='TEXT'; text=put(ix,2.-L); style="SWISS"; size=1; x=x2+0.25; output;
*/
  format stage 4.;
  keep x y xsys ysys function stage color size ;*text style ;
run;

%let pct = %sysevalf (100/(&nSets+.5));
symbol1  h=&pct.pct v=plus    color=CX505050;
symbol2  h=&pct.pct v=square  color=blue;
symbol3  h=&pct.pct v=circle  color=green;
symbol4  h=&pct.pct v=diamond color=red;
symbol5  h=&pct.pct v=square  color=green;
symbol6  h=&pct.pct v=circle  color=red;
symbol7  h=&pct.pct v=diamond color=blue;
symbol8  h=&pct.pct v=square  color=red;
symbol9  h=&pct.pct v=circle  color=blue;
symbol10 h=&pct.pct v=diamond color=green;
axis1 order=(0 to %eval(1+&nSets)) minor=none label=none value=none major=none;
axis2 order=(0 to 27)    minor=none  label=none value=none major=none;
legend1 label=none value=none;

*goptions gsfname=gout gsfmode=replace device=gifanim delay=60;
*goptions hsize=8in vsize=6in;
*goptions hsize=4in vsize=3in;
*filename gout "c:\temp\gif\superset.gif";
*ods listing ;
*uncomment these to create individual gifs named stage001-stageNNN;
*the individual gifs can be into a supercompressed gif animation
*using Alchemy Mindworks GIF Construction Set;
goptions device=gif;
goptions hsize=8.42in vsize=6.31in; * 800x600;
goptions hsize=4.21in vsize=3.15in; * 400x300;
ods listing close;
ods html gpath="c:\temp\gif" body="super.html";
goptions i=none goutmode=replace ftext='Arial' htext=5pct hby=0;
*title h=5pct move=(1pct,1pct) "#byvar1 #byval1";
data stages;
  set stages;
  y = &nSets - y + 1;
run;
data annotate;
  set annotate;
  y = &nSets - y + 1;
run;
proc gplot data=stages uniform ;
  by stage;
  plot y * x = state
  / vaxis=axis1
    haxis=axis2
/*
    autohref lautohref=33 cautohref=cxAAAAAA
    autovref lautovref=33 cautovref=cxAAAAAA
*/
    nolegend
    anno=annotate
    name="stage001"
    noframe noaxis
  ;
*  where stage<30;
run;
quit;
ods html close;
data _null_;
   file gout recfm=n mod;
   put '3B'x;
run;
goptions reset=all;

options noxwait noxsync;
x "explorer c:\temp\gif\superset.gif";
*/;
options source;