options ps=500 nodate;

title '   ';
PROC IML;

print 'BLOCKING USING DS';

start DESIGN;



   TRIES=100;      /* Enter Number of Tries */
   BLOCK=4;        /* Enter Total Number of Blocks */
   NB=3;           /* Enter Number of High Score Blocks (a) */
   BSIZE={9 9 9 9};/* Enter Vector of Block 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 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 1 1 1 1 1};

   print TRIES BLOCK NB ;
   print WEI;

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

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

/* B is the design matrix for blocks */

   B=j(BSIZE[1],1,1);
   do I=2 to BLOCK;
      B=B//j(BSIZE[I],1,I);
   end;
   B=design(B);

finish;

start SWAP;

  BET2=0;
  do J=1 to NB;
     BL1=BLOC[J];
     BS1=BSIZE[BL1];
     AB=0;
     do I=1 to NB;
        if AB>0 | I>J then do;
           BL2=BLOC[I];
           BS2=BSIZE[BL2];
           do K1=KI[BL1,1] to KI[BL1,2];
              do K2=KI[BL2,1] to KI[BL2,2];
                 if any(X[K1,2:KFAC1]=X[K2,2:KFAC1]) then do;
                    XP=X#B1[,BL2];
                    XQ=X#B1[,BL1];
                    A2=(X[K2,]-X[K1,])`;
                    A1=(XP[+,]/BS2-XQ[+,]/BS1)`-A2;
                    BKVU=(A2||(A1+((BS1-1)/BS1+(BS2-1)/BS2)#A2))`*XINVU;
                    IT=I(2)+BKVU*(A1||A2);
                    if det(IT)^=0 then do;
                       VNC=XINVU-XINVU*(A1||A2)*inv(IT)*BKVU;
                       CRIN=det(VNC[2:P1,2:P1]);
                       if CRIN<CRITERI then do;
                          CRITERI=CRIN;
                          XI=X;
                          XI[(K2//K1),]=X[(K1//K2),];
                          X=XI;
                          BET=1;
                          XINVU=VNC;
                          AI=XINV*A2;
                          SCORB0[,BL1]=SCORB0[,BL1]+AI;
                          SCORB0[,BL2]=SCORB0[,BL2]-AI;
                          BET2=1;
                          AB=1;
                       end;
                    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`*PP*X;
      DETI=det(IT);

finish;

start SCORE;

/* FINDS THE PARAMETERS THAT HAVE THE SMALLEST EFFICIENCIES */
/* AND FINDS THE BLOCKS THAT HAVE THE HIGHEST SCORES FOR THAT PARAMETERS */

      BET=0;
      MSCOR=abs(SCORB[+,]);
      BLOC=shape(0,1,NB);
      MSCOR=rank(MSCOR);
      do JJ=1 to NB;
         BLOC[JJ]=loc(MSCOR=BLOCK1-JJ);
      end;
      BET2=1;
      do until(BET2=0);
         run SWAP;
      end;



finish;

start SEARCH;

/* CONTROLS THE SEARCH PROCEDURE */

   CRIBE=100;
   TENTR=0;
   do TRY=1 to TRIES;
      DETI=0;
      do until(DETI>0);
         run GENERATE;
      end;
      XINVU=inv(IT);
      XINVU=XINVU[1:P1,1:P1];
      CRITERI=det(XINVU[2:P1,2:P1]);
      SCOR=(XINV*X[1:N,]`);
      SCORB0=SCOR*B;
      BET=1;
      do until(BET=0);
         SCORB=SCORB0[2:P1,];
         run SCORE;
      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;
      SCORBL=SCORB0;

finish;

start RESULT;

/* CALCULATES SOME QUANTITIES FOR THE FINAL DESIGN */
/* AND PRINTS THE RESULTS */
      PRINT FTRY;
      if TENTR=1 then do;
         DC=DATAC;
         DATAC[rank(DATAC)]=DC;
         DATAC=DATAC[1:10];
      end;
      BEST=BESTX[,2:KFAC1];
      SCORB=BESTX`*B-(1/N)*BESTX`*j(N,N,1)*B;
      SSS=SCORBL#SCORBL;
      SSS=(SSS[,+]);
      SSS=WEI*SSS[2:P1]/sum(WEI);
      SCORB=SCORB[2:P1,];
      NYS=SCORB#SCORB;
      NYS=(NYS[,+]);
      AS=WEI*VARB/sum(WEI);
      EFVB=VAR/VARB;
      WMEF=WEI*EFVB/sum(WEI);
      NYG=WEI*NYS/sum(WEI);
      print DATAC;
      print 'RESULTING DESIGN';
      print TRIES BLOCK NB;
      print WEI;
      print BEST;
      print CRIBE SSS WMEF AS NYG;
      print VARB EFVB;

finish;



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

  KI=shape(0,BLOCK,2);
  A=0;
  do I=1 to BLOCK;
     KI[I,1]=A+1;
     KI[I,2]=A+BSIZE[I];
     A=KI[I,2];
  end;

/* 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(1,P1,0);
  do L3=1 to BLOCK;
     CON=CON||BSIZE[L3]/N;
  end;
  B1=B//CON;
  INVB=inv(B1`*B1);
  PP=i(N+1)-B1*INVB*B1`;
  run SEARCH;
  run RESULT;

quit;