options ps=500 nodate;
title '   ';
PROC IML;

/* This Program is for Selecting the Best Design Taking into Account
   the Interactions Between Lower and Higher Stratum Factors */


print 'INTERACTION PHASE OF MULTILEVEL STRUCTURE DESIGN';


start DESIGN;

   TRIES=10;         /* Enter Number of Tries       */
   BLOCK=6;         /* Enter Total Number of Blocks */
   NF1=1;            /* Enter the Number of Factors in the Previous Level*/
   BSIZE={3,3,3,3,3,3};          /* Enter Vector of Block Sizes */

/* Enter Presence (1)/ Absence (0)  */
/* of Terms in the Model - Based on */
/* the Full Second-order Model      */
/* Main Effects of the Higher Stratum Factors should not be in the Model */

   MODEL={0 1 1 0 1 1 1 1 1};


/* Enter Weights for each Interaction.          */

   WEI={1 1 1};

/* Enter the Treatment Set */
/* XX is the treatment set Matrix */

XX={ 0  0  0  1,
     0  0  1  0,
    -1  0  0  0,
     0  0  0  0,
     0  1  0  0,
     0  0  0  0,
    -1  1 -1  1,
     1  1  1 -1,
     1 -1 -1 -1,
    -1 -1  1  1,
     0  0  0  0,
     1  1  1  1,
    -1  1 -1 -1,
    -1 -1  1 -1,
     1 -1 -1  1,
     0  0 -1  0,
     1  0  0  0,
     0  0  0 -1,
     0 -1  0  0,
     0  0  0  0,
     0  0  0  0,
     1  1 -1  1,
     1 -1  1 -1,
    -1  1  1 -1,
    -1 -1 -1  1,
     1 -1  1  1,
     0  0  0  0,
     1  1 -1 -1,
    -1  1  1  1,
    -1 -1 -1 -1};

XN=XX;

/* B is the design matrix for blocks */

   B=i(BLOCK)@j(BSIZE,1,1);

finish;

start SWAP;


  BET=0;
  do J=1 to BLOCK;
     AB=0;
     R1=(J-1)*BSIZE;
     do I=1 to BLOCK;
        if AB>0 | I>J then do;
           DIF=0;
           R2=(I-1)*BSIZE;
           if any(XN[R1+1:R1+BSIZE,1:NF1]^=XN[R2+1:R2+BSIZE,1:NF1]) then do;
              XI=XN;
              XI[R1+1:R1+BSIZE,NF1+1:KFAC]=XN[R2+1:R2+BSIZE,NF1+1:KFAC];
              XI[R2+1:R2+BSIZE,NF1+1:KFAC]=XN[R1+1:R1+BSIZE,NF1+1:KFAC];
              run MODEL(XI);
              IT=XM`*PP*XM;
              if det(IT)^=0 then do;
                 VNC=inv(IT);
                 VNC=vecdiag(VNC);
                 CRIN=WEI*VNC/sum(WEI);
                 if CRIN<CRITERI then do;
                    CRITERI=CRIN;
                    XN=XI;
                    X=XI;
                    VARF=VNC;
                    BET=1;
                    AB=1;
                 end;
              end;
           end;
        end;
     end;
  end;

finish;

start MODEL(XA) global(KFAC,XX,NF1,XM,MODEL);

      free XM;
      XM=XA||XA#XA;
      do II=1 to NF1;
         do JJ=II+1 to KFAC;
            XM=XM||XX[,II]#XA[,JJ];
         end;
      end;
      do II=NF1+1 to KFAC-1;
         do JJ=II+1 to KFAC;
            XM=XM||XA[,II]#XA[,JJ];
         end;
      end;
      XM=XM[,loc(MODEL)];

finish;

start GENERATE;

/* SELECTS A RANDOM INITIAL DESIGN */

      free XN;
      Y=uniform(j(BLOCK,1,0));
      YR=rank(Y);
      do I=1 to BLOCK;
         R=(YR[I]-1)*BSIZE;
         XN=XN//XX[R+1:R+BSIZE,NF1+1:KFAC];
      end;
      XN=XX[,1:NF1]||XN;
      run MODEL(XN);
      IT=XM`*PP*XM;
      DETI=det(IT);

finish;



start SEARCH;

/* CONTROLS THE SEARCH PROCEDURE */

   CRIBE=10000;
   do TRY=1 to TRIES;
      DETI=0;
      do until(DETI>0);
         run GENERATE;
      end;
      XINVU=inv(IT);
      VARF=vecdiag(XINVU);
      CRITERI=WEI*VARF/sum(WEI);
      BET=1;
      do until(BET=0);
         run SWAP;
      end;
      if CRITERI<CRIBE then run FINTRY;
   end;
   run RESULT;

finish;



start FINTRY;

/* KEEPS THE RESULTS FOR ONE TRY */

      BEST=X;
      CRIBE=CRITERI;
      VARB=VARF;

finish;

start RESULT;

/* CALCULATES SOME QUANTITIES FOR THE FINAL DESIGN */
/* AND PRINTS THE RESULTS */

      print TRIES  BLOCK;
      print WEI;
      print BEST;
      print CRIBE;
      print VARB ;
      stop;

finish;

  reset storage='SPLIT';
  load BEST;

  run DESIGN;
  KFAC=ncol(XX);
  P=sum(MODEL);
  N=nrow(B);
  INVB=inv(B`*B);
  PP=i(N)-B*INVB*B`;
  run SEARCH;
  reset storage='SPLIT2';
  store BEST;

quit;