options ps=500 nodate;

title '   ';
PROC IML;

start DESIGN;

   TRIES=1;        /* Enter Number of Tries */
   ROW=7;            /* Enter Total Number of Rows */
   COL=4;            /* Enter Total Number of Columns */

/* 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 Linear, */
/* for each Quadratic,            */
/* for each Interaction.          */

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

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

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


/* R and C are the design matrices for Rows and Columns */

   R=i(ROW)@j(COL,1,1)//j(1,ROW,1)//j(1,ROW,0);
   C=j(ROW,1,1)@i(COL)//j(1,COL,0)//j(1,COL,1);
   B=C||R;
   PP=B*inv(B`*B)*B`;

finish;

start COLUMN;

   STOPC=1;
   if JC^=IC then do;
      XP=X#C[,JC];
      XQ=X#C[,IC];
      A1=(XP[+,]-XQ[+,])`/ROW-A2;
      AC=A1||A2;
      BKVC=(A2||(A1+((ROW-1)/ROW+(ROW-1)/ROW)#A2))`*VNR;
      IT=I(2)+BKVC*AC;
      if det(IT)>0 then do;
         VNC=VNR-VNR*AC*inv(IT)*BKVC;
         STOPC=0;
      end;
   end;
   else do;
      VNC=VNR;
      STOPC=0;
   end;

finish;

start SWAP;

  BET=0;
  do J=1 to N;
     do I=1 to N;
        if any(X[J,]^=X[I,]) then do;
           JR=int(J/COL-1/(COL*2))+1;
           IR=int(I/COL-1/(COL*2))+1;
           JC=mod(J-1,COL)+1;
           IC=mod(I-1,COL)+1;
           A2=(X[J,]-X[I,])`;
           STOPR=1;
           if JR^=IR then do;
              XP=X#R[,JR];
              XQ=X#R[,IR];
              A1=(XP[+,]-XQ[+,])`/COL-A2;
              AR=A1||A2;
              BKVR=(A2||(A1+((COL-1)/COL+(COL-1)/COL)#A2))`*XINVU;
              IT=I(2)+BKVR*AR;
              if det(IT)>0 then do;
                 VNR=XINVU-XINVU*AR*inv(IT)*BKVR;
                 STOPR=0;
              end;
           end;
           else do;
              VNR=XINVU;
              STOPR=0;
           end;
           if STOPR=0 then do;
              run COLUMN;
              if STOPC=0 then do;
                 VNF=vecdiag(VNC[2:P1,2:P1]);
                 EFVARN=(VAR/VNF)*100;
                 CRIN=WEI*EFVARN/sum(WEI);
                 if CRIN>CRITERI then do;
                    CRITERI=CRIN;
                    XI=X;
                    XI[(I//J),]=X[(J//I),];
                    IT2=(r||c||xi)`*(r||c||xi);
                    if DET(IT2)>0 then do;
                       XINVU=VNC;
                       X=XI;
                       EFFVARF=EFVARN;
                       BET=1;
                    end;
                 end;
              end;
           end;
        end;
     end;
  end;
finish;

start MODEL;

      XT=XX||XX#XX;
      do I=1 to KFAC-1;
         do J=I+1 to KFAC;
            XT=XT||XX[,I]#XX[,J];
         end;
      end;
      X=j(N,1)||XT[,loc(MODEL)];

finish;

start GENERATE;

/* SELECTS A RANDOM INITIAL DESIGN */

      Y=uniform(j(N,1,0));
      YR=rank(Y);
      X=X[YR,];
      X=X//CON1;
      IT=X`*(i(N+2)-PP)*X;
      DETI=det(IT);

finish;

start SEARCH;

/* CONTROLS THE SEARCH PROCEDURE */

   CRIBE=0;
   TENTR=0;
   do TRY=1 to TRIES while(CRIBE<99.999999);
      DETI=0;
      do until(DETI>0);
         run GENERATE;
      end;
      XINVU=inv(IT);
      VARI=vecdiag(XINVU[2:P1,2:P1]);
      EFFVARF=(VAR/VARI)*100;
      CRITERI=WEI*EFFVARF/sum(WEI);
      BET=1;
      do until(BET=0);
         run SWAP;
      end;
      DATAC=DATAC//CRITERI;
      if CRITERI>CRIBE then run FINTRY;
      if TRY=10 then TENTR=1;
      FTRY=TRY;
   end;

finish;


start FINTRY;

/* KEEPS THE RESULTS FOR ONE TRY */

      BESTX=X[1:N,];
      XINVF=XINVU;
      VARB=vecdiag(XINVF[2:P1,2:P1]);
      CRIBE=CRITERI;
      EFVB=EFFVARF;

finish;

start RESULT;

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

      if TENTR=1 then do;
         DC=DATAC;
         DATAC[FTRY+1-RANK(DATAC)]=DC;
         DATAC=DATAC[1:10];
      end;
      print DATAC;
      C=C[1:N,]//j(1,COL,1);
      R=R[1:N,]//j(1,ROW,1);
      PP=i(N+1)-C*inv(C`*C)*C`;
      X=BESTX//j(1,P1,0);
      XINC=inv(X`*PP*X);
      VARC=vecdiag(XINC[2:P1,2:P1]);
      EFFC=(VAR/VARC)*100;
      CRITC=WEI*EFFC/sum(WEI);
      PP=i(N+1)-R*inv(R`*R)*R`;
      XINC=inv(X`*PP*X);
      VARR=vecdiag(XINC[2:P1,2:P1]);
      EFFR=(VAR/VARR)*100;
      CRITR=WEI*EFFR/sum(WEI);
      BEST=BESTX[,2:KFAC1];
      AS=WEI*VARB/sum(WEI);
      DS=det(XINVF[2:P1,2:P1]);
      print 'RESULTING DESIGN';
      print TRIES ROW COL;
      print WEI;
      print BEST;
      print 'Column Design';
      print COL CRITC;
      print VARC EFFC;
      print 'Row Design';
      print ROW CRITR;
      print VARR EFFR;
      print 'Combined Design';
      print ROW COL CRIBE DS AS;
      print VARB EFVB;

finish;



  run DESIGN;
  KFAC=ncol(XX);
  KFAC1=KFAC+1;
  P1=sum(MODEL)+1;
  N=nrow(XX);
  Y=j(N,1,0);

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

  run MODEL;
  XINV=INV(X`*X);
  VAR=vecdiag(XINV[2:P1,2:P1]);
  CON1=j(2,P1,0);
  run SEARCH;
  run RESULT;

quit;