/* 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 [i] = 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 [i] = 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;