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

print 'LAST PHASE OF MULTILEVEL STRUCTURE DESIGN';

   reset storage='SPLIT';
   load BEST;

start DESIGN;

   TRIES=10;        /* Enter Number of Tries       */
   BLOCK=2;         /* Enter Total Number of Blocks */
   NF1=0;           /* Enter the Number of Factors at the Previous Stratum */

   BSIZE={3,3};              /* Enter Vector of Block Sizes */
   MAIN=6;              /* Enter Total Number of Main Units */
   MSIZE={3,3,3,3,3,3};             /* Enter Vector of Main Unit Sizes */

/* Enter Presence (1)/ Absence (0)  */
/* of Terms in the Model - Based on */
/* the Full Second-order Model      */
/* The Intercept is included automatically*/

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


/* Enter Weights for each Interaction.          */

   WEI={1 1 1 1 1 1 1 1 1};

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

XX={
};
X=XX;

/* B is the design matrix for blocks */

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

finish;

start SWAP;


  BET=0;
  do J=1 to MAIN;
     AB=0;
     R1=(J-1)*MSIZE;
     do I=1 to MAIN;
        if AB>0 | I>J then do;
           R2=(I-1)*MSIZE;
           AX=j(MSIZE,KFAC,0);
           BX=j(MSIZE,KFAC,0);
           if NF1>0 then do;
              AX=X[R1+1:R1+MSIZE,1:NF1];
              BX=X[R2+1:R2+MSIZE,1:NF1];
           end;
           if AX=BX | NF1=0 then do;
              XI=X;
              XI[R1+1:R1+MSIZE,]=X[R2+1:R2+MSIZE,];
              XI[R2+1:R2+MSIZE,]=X[R1+1:R1+MSIZE,];
              run MODEL(XI);
              XM=XM//CON1;
              IT=XM`*PP*XM;
              if det(IT)^=0 then do;
                 VNC=inv(IT);
                 VNC=vecdiag(VNC[2:P1,2:P1]);
                 EFFVARN=(VAR/VNC)#100;
                 CRIN=WEI*EFFVARN/sum(WEI);
                 if CRIN>CRITERI then do;
                    CRITERI=CRIN;
                    EFFVARF=EFFVARN;
                    X=XI;
                    VARF=VNC;
                    BET=1;
                    AB=1;
                 end;
              end;
           end;
        end;
     end;
  end;

finish;

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

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

finish;

start GENERATE;

/* SELECTS A RANDOM INITIAL DESIGN */

      XN=XX[,NF1+1:KFAC];
      do I=1 to MAIN-1;
         do J=2 to MAIN;
            R1=(I-1)*MSIZE;
            R2=(J-1)*MSIZE;
            AX=j(MSIZE,KFAC,0);
            BX=j(MSIZE,KFAC,0);
            if NF1>0 then do;
               AX=XX[R1+1:R1+MSIZE,1:NF1];
               BX=XX[R2+1:R2+MSIZE,1:NF1];
            end;
            if AX=BX | NF1=0 then do;
               U=uniform(0);
               if U>0.5 then do;
                  XN[R1+1:R1+MSIZE,]=XX[R2+1:R2+MSIZE,NF1+1:KFAC];
                  XN[R2+1:R2+MSIZE,]=XX[R1+1:R1+MSIZE,NF1+1:KFAC];
               end;
            end;
         end;
      end;
      if NF1>0 then XN=XX[,1:NF1]||XN;
      run MODEL(XN);
      XM=XM//CON1;
      IT=XM`*PP*XM;
      DETI=abs(det(IT));

finish;



start SEARCH;

/* CONTROLS THE SEARCH PROCEDURE */

   CRIBE=0;
   do TRY=1 to TRIES while(CRIBE<99.999999);
      DETI=0;
      do until(DETI>0);
         run GENERATE;
      end;
      XINVU=inv(IT);
      VARF=vecdiag(XINVU[2:P1,2:P1]);
      EFFVARF=(VAR/VARF)#100;
      CRITERI=WEI*EFFVARF/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;
      EFVB=EFFVARF;

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;

  run DESIGN;
  KFAC=ncol(XX);
  P=sum(MODEL);
  N=nrow(B);
  P1=P+1;


/* GENERATES THE COMPLETE X MATRIX AND CALC. VAR BEFORE BLOCKING */

  run MODEL(XX);
  XINV=INV(XM`*XM);
  VAR=vecdiag(XINV[2:P1,2:P1]);
  CON1=j(1,P1,0);
  CON=j(1,BLOCK,BSIZE/N);
  B1=B//CON;
  PP=i(N+1)-B1*inv(B1`*B1)*B1`;
  run SEARCH;
  reset storage='SPLIT3';
  store BEST;

quit;