/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 21.3                                              */
/*   TITLE: Example 21.2 and SAS Task 21.2-21.3: Response surface in  */
/*          four factors                                              */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Uses an IML-implemented exchange search.  Took ~9 minutes */
/*          to run on server-class PC in May, 2007.                   */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 21.2: Response surface in four factors";

data Can;
   do x1 = -1 to 1 by 1;
   do x2 = -1 to 1 by 1;
   do x3 = -1 to 1 by 1;
   do x4 = -1 to 1 by 1;
      output;
      end; end; end; end;
run;

proc iml;
   use Can;
   read all var {x1 x2 x3 x4} into xCan;

   start MakeZ(x);
      F1 = j(nrow(x),1) || x
           || x[,1]#x[,2] || x[,1]#x[,3] || x[,1]#x[,4]
                          || x[,2]#x[,3] || x[,2]#x[,4]
                                         || x[,3]#x[,4];
      F2 = x[,1]#x[,1] || x[,2]#x[,2] || x[,3]#x[,3] || x[,4]#x[,4];
      F = F1 || F2;
      return(F);
   finish;

   eps = 1e-8;
   start DCritW(w) global(F,k,eps);
      p = ncol(F);
      r = 11;
      s = p - r;
      F1 = F[,1:r];
      Dw = diag((w##2)/sum(w##2));
      m11 = F1`*Dw*F1 + eps*i(ncol(F1));
      M   = F `*Dw*F  + eps*i(ncol(F ));
      return((k/r)*log(det(M11)) + ((1-k)/s)*(log(det(M)) - log(det(M11))));
   finish;

   start KEval(w) global(xCan);
      free w0 w1 w2 w3 w4;
      do i = 1 to nrow(xCan);
         nnz = sum(xCan[i,] ^= 0);
         if      (nnz = 0) then w0 = w0 // w[i];
         else if (nnz = 1) then w1 = w1 // w[i];
         else if (nnz = 2) then w2 = w2 // w[i];
         else if (nnz = 3) then w3 = w3 // w[i];
         else                   w4 = w4 // w[i];
         end;
      return(   (w0[+] || sqrt(((w0 - w0[:])##2)[:]))
             // (w1[+] || sqrt(((w1 - w1[:])##2)[:]))
             // (w2[+] || sqrt(((w2 - w2[:])##2)[:]))
             // (w3[+] || sqrt(((w3 - w3[:])##2)[:]))
             // (w4[+] || sqrt(((w4 - w4[:])##2)[:])));
   finish;

   /*
   /  Ds-optimum design
   /------------------------------------------------------------------*/
   k = 0;
   x     = XCan; F = MakeZ(x);
   winit = j(nrow(X),1);
   call nlpqn(rc,wOp,"DCritW",winit,1);
   wOp = ((wOp##2)/sum(wOp##2))`;
   XOp0 = X  [loc(wOp > 1e-3),];
   wOp0 = wOp[loc(wOp > 1e-3),];

   /*
   /  D-optimum design
   /------------------------------------------------------------------*/
   k = 1;
   x     = XCan; F = MakeZ(x);
   winit = sqrt(j(nrow(X),1));
   call nlpqn(rc,wOp,"DCritW",winit,1);
   wOp = ((wOp##2)/sum(wOp##2))`;
   XOp1 = X  [loc(wOp > 1e-3),];
   wOp1 = wOp[loc(wOp > 1e-3),];

   start D1(w,x) global(wOp0,xOp0);
      F = MakeZ(x);
      p = ncol(F);
      r = 11;
      s = p - r;
      F1 = F[,1:r];
      Dw = diag(w);
      m11 = F1`*Dw*F1;
      if (det(M11)<=0) then
         d1 = -999;
      else
         d1 = log(det(m11))/r;
      return(d1);
   finish;

   start Ds(w,x) global(wOp1,xOp1);
      F = MakeZ(x);
      p = ncol(F);
      r = 11;
      s = p - r;
      F1 = F[,1:r];
      Dw = diag(w);
      m11 = F1`*Dw*F1;
      M   = F `*Dw*F ;
      if (det(M) <= 0) then
         ds = -999;
      else if (det(M11) <= 0) then
         ds = 999;
      else
         ds = (log(det(M)) - log(det(M11)))/s;
      return(ds);
   finish;

   start EffD(w,x) global(wOp1,xOp1);
      d   = D1(w   ,x   );
      dOp = D1(wOp1,xOp1);
      return(exp(d)/exp(dOp));
   finish;

   start EffDs(w,x) global(wOp0,xOp0);
      d   = Ds(w   ,x   );
      dOp = Ds(wOp0,xOp0);
      return(exp(d)/exp(dOp));
   finish;

   /*
   /  D/Ds-optimum continuous design
   /------------------------------------------------------------------*/
   options nonotes;
   k = 0.5;
   x = XCan;
   F = MakeZ(x);
   winit = sqrt(j(nrow(X),1));
   call nlpqn(rc,w,"DCritW",winit,1);
   wOp = w;
   wOp = ((wOp##2)/sum(wOp##2))`;
   XOpC = X  [loc(wOp > 1e-3),];
   wOpC = wOp[loc(wOp > 1e-3),];

   F = MakeZ(xOpC);
   dCOp = DCritW(sqrt(wOpC));

   start DCritX(x) global(k);
      F = MakeZ(shape(x,nrow(x)*ncol(x)/4,4));
      p = ncol(F);
      r = 11;
      s = p - r;
      F1 = F[,1:r];
      m11 = F1`*F1/nrow(F1) + 1e-8*i(ncol(F1));
      M   = F `*F /nrow(F ) + 1e-8*i(ncol(F ));
      return((k/r)*log(det(M11)) + ((1-k)/s)*(log(det(M)) - log(det(M11))));
   finish;


   /*
   /  Exchange search for best design on grid.
   /------------------------------------------------------------------*/
   if (inputn(symget('sysver'),'3.1') >= 9) then do;
      start sort(x,idx);
         sortx = j(nrow(x), ncol(x), .);
         do i = nrow(idx)*ncol(idx) to 1 by -1;
            sortx[rank(x[,idx[i]]),] = x[];
            end;
         x = sortx;
      finish;
      end;
   start Exchange;
      free AllDesigns;
      bestDCritX = -1e10;
      do i = 1 to 1000;
         xCurr =    int(3*ranuni(1:N)`) - 1 
                 || int(3*ranuni(1:N)`) - 1
                 || int(3*ranuni(1:N)`) - 1
                 || int(3*ranuni(1:N)`) - 1;

         dCurr = dCritX(xCurr);
         lastD = dCurr - 1;
         iAddDel = 0;
         do while (dCurr > lastD+1e-8);
            lastD = dCurr;

            bestdAdd = -1e10;
            do iAdd = 1 to nrow(xCan);
               xAdd = xCurr // xCan[iAdd,];
               dAdd = dCritX(xAdd);
               if (dAdd > bestdAdd) then do;
                  bestdAdd = dAdd;
                  bestiAdd = iAdd;
                  end;
               end;
            xCurr = xCurr // xCan[bestiAdd,];

            bestdDel = -1e10;
            do iDel = 1 to nrow(xCurr);
               xDel = xCurr[loc((1:nRow(xCurr)) ^= iDel),];
               dDel = dCritX(xDel);
               if (dDel > bestdDel) then do;
                  bestdDel = dDel;
                  bestiDel = iDel;
                  end;
               end;
            xCurr = xCurr[loc((1:nRow(xCurr)) ^= bestiDel),];

            dCurr = dCritX(xCurr);
            iAddDel = iAddDel + 1;
            end;
         xOp = xCurr;
         if (DCritX(xOp) > bestDCritX) then do;
            bestDCritX = DCritX(xOp);
            bestxOp    = xOp;
            end;
         AllDesigns = AllDesigns // (i || DCritX(xOp) || bestDCritX);
         end;
      xOp = bestxOp;
      dOp = DCritX(xOp);
      AllDesigns = -AllDesigns;
      call sort(AllDesigns,2);
      AllDesigns = -AllDesigns;

/*
/  Uncomment to see how well we're doing with the optimization.
/---------------------------------------------------------------------
      print (AllDesigns[1:10,]);
*/

      i = ((xOp+1)*(3##(3:0))`)+1;
      r = rank(i);
      j = i;
      i[r] = j;
      xOpG = xCan[i,];
      dOpG = DCritX(xOpG);
      wOpG = shape(1/nrow(xOpG),nrow(xOpG),1);

      /*
      /  Nonlinear optimization for best design off-grid.
      /---------------------------------------------------------------*/
      xinit = shape(xOpG,1,nrow(xOpG)*ncol(xOpG));
      con =    (shape(-1,1,N) || shape(-1,1,N) || shape(-1,1,N) || shape(-1,1,N))
            // (shape( 1,1,N) || shape( 1,1,N) || shape( 1,1,N) || shape( 1,1,N));
      call nlpqn(rc,x,"DCritX",xinit,1,con);
      xOpO = shape(x,nrow(x)*ncol(x)/4,4);
      dOpO = DCritX(xOpO);
      wOpO = shape(1/nrow(xOpO),nrow(xOpO),1);

      Compound = (exp(dOpG)        // exp(dOpO)       ) / exp(dCOp);
	  EffD     =  EffD (wOpG,xOpG) // EffD (wOpO,xOpO);
	  EffDs    =  EffDs(wOpG,xOpG) // EffDs(wOpO,xOpO);
      print N Compound [format=5.3] EffD [format=5.3] EffDs [format=5.3];
   finish;

   print "Table 21.3 - Efficiencies for exact compound optimum designs";

   N = 15; call Exchange;
   N = 20; call Exchange;
   N = 25; call Exchange;

   Compound = exp(dCOp) / exp(dCOp);
   EffD     = EffD (wOpC,xOpC);
   EffDs    = EffDs(wOpC,xOpC);
   print (.) Compound [format=5.3] EffD [format=5.3] EffDs [format=5.3];
quit;

title1;