/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 19.3                                              */
/*   TITLE: Tables 19.2 & 19.3 - Second order response surface,       */
/*          constrained design region                                 */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Second order response surface, constrained design region";

/*
/  Create vertices of the design region.
/---------------------------------------------------------------------*/
data Candidates;
   do x1 = -1 to 1 by 0.1;
      do x2 = -1 to 1 by 0.1;
         output;
         end;
      end;
data Candidates; set Candidates;
   x1 = round(x1,0.1);
   x2 = round(x2,0.1);
   if (  (2*x1 + x2 <= 1       + 1e-6)
       & (-1        <= x1 + x2 + 1e-6)
       & (  x2 - x1 <= 1.5     + 1e-6)) then output;
run;

title2 "First-order D-optimum design";
proc optex data=Candidates;
   model x1 x2;
   generate n=3 method=m_fedorov niter=1000 keep=10;
   output out=FirstOrder;
proc print data=FirstOrder;
run;
title2;

/*
/  Compute optimum weights for second-order model over candidates.
/---------------------------------------------------------------------*/
proc iml;
   start MakeZ(_x);
      x = shape(_x,ncol(_x)*nrow(_x)/2,2);
      Z = j(nrow(x),1) || x[,1] || x[,2]
                       || x[,1]#x[,1]
                       || x[,1]#x[,2]
                       || x[,2]#x[,2];
      return(Z);
   finish;

   start llog(x); if (x > 0) then return(log(x)); else return(-9999); finish;

   start DCritW(w) global(Z);
      d = det(Z`*diag((w##2)/sum(w##2))*Z);
      return(llog(d)/ncol(Z));
   finish;

   use Candidates; read all var {x1 x2} into xCan;

   x     = xCan; Z = MakeZ(x);
   winit = sqrt(j(nrow(x),1) / nrow(x));
   call nlpqn(rc,w,"DCritW",winit,1);
   w = ((w##2)/sum(w##2))`;
   xOp = x[loc(w > 1e-5),];
   wOp = w[loc(w > 1e-5) ]; wOp = wOp /sum(wOp);

   x     = xOp; Z = MakeZ(x);
   winit = sqrt(wOp);
   call nlpqn(rc,w,"DCritW",winit,1);
   w = ((w##2)/sum(w##2))`;
   xOp = x[loc(w > 1e-5),];
   wOp = w[loc(w > 1e-5) ]; wOp = wOp /sum(wOp);

   wEx = round(19*wOp); wEx = wEx / sum(wEx);
   Z = MakeZ(xOp);
   DEffExact =   exp(DCritW(sqrt(wEx)))
               / exp(DCritW(sqrt(wOp)));
   call symput('DEffExact',trim(left(putn(DEffExact,"percent9.2"))));

   data = wOp || xOp;
   create SecondOrder var {w x1 x2}; append from data;
run;

data FirstOrder; set FirstOrder; w = 1 / 3;
proc sort data=FirstOrder;  by x2 x1;
proc sort data=SecondOrder; by x2 x1;
data Table1902;
   merge FirstOrder (rename=(w=w1))
         SecondOrder(rename=(w=w2));
   by x2 x1;
   r = round(19*w2);
   label w1 = "First-Order Weight"
         w2 = "Second-Order Weight"
         r  = "[19*w]";
run;

title2 "Table 19.2 - D-optimum first- and second-order designs";
title3 "D-efficiency of exact design = &DEffExact";
proc print data=Table1902 noobs label;
   format w1 w2 6.4;
run;
title3;
title2;

/*
/  Compute 9+1 design optimum for checking the 9-run D-optimum
/  first-order design.
/---------------------------------------------------------------------*/
title2 "9-run D-optimum first-order design";
proc optex data=Candidates coding=orthcan;
   model x1 x2;
   generate n=9 method=m_fedorov niter=1000 keep=10;
   output out=D1;
run;
title2 "9+1 design optimum for augmenting first-order design to check model";
proc optex data=Candidates coding=orthcan;
   model x1 x2,x1|x1|x2|x2@2 / prior=0,1;
   generate n=10 method=m_fedorov augment=D1 niter=1000 keep=10;
   output out=D2;
proc freq data=D2 noprint;
   table x2*x1 / list out=Check(rename=(Count=rCheck));
run;
title2;


/*
/  Select points for augmenting the first-order design from the
/  optimum second-order design.
/---------------------------------------------------------------------*/
data Initial;    set Check;       Set = "Initial   ";
   do i = 1 to rCheck; output; end;
data Candidates; set SecondOrder; Set = "Candidates";
run;

proc iml;
   use Initial;    read all var {x1 x2};
   n0    = nrow(x1);
   FInit = j(n0,1) || x1 || x2 || x1#x1 || x1#x2 || x2#x2;

   use Candidates; read all var {x1 x2};
   nCan  = nrow(x1);
   FCan  = j(nCan,1) || x1 || x2 || x1#x1 || x1#x2 || x2#x2;

   start DCritW(w) global(FInit,F,alpha,n0,nCan);
      Dw = diag((w##2)/sum(w##2));
      return(log(det(  (1-alpha)*FInit`   *FInit/n0
                     +    alpha *F    `*Dw*F       )));
   finish;

   start AlphaOpt;
      F     = FCan;
      winit = sqrt(j(nCan,1) / nCan);
      call nlpqn(rc,w,"DCritW",winit,1);
      w = ((w##2)/sum(w##2))`;

      F     = F[loc(w > 1e-5),];
      winit = sqrt(w[loc(w > 1e-5),]);
      call nlpqn(rc,w,"DCritW",winit,1);
      w = ((w##2)/sum(w##2))`;

      FOp   = F[loc(w > 1e-5),];
      wOp   = w[loc(w > 1e-5),];
      wOp   = wOp/sum(wOp);
   finish;

   free FAll wAll Which;

   alpha = 5/(10 + 5);
   call AlphaOpt;
   data = FOp[,2:3] || wOp;
   create Add333 var {x1 x2 w};
   append from data;

   alpha = 9/(10 + 9);
   call AlphaOpt;
   data = FOp[,2:3] || wOp;
   create Add474 var {x1 x2 w};
   append from data;
quit;

proc sort data=Check ; by x2 x1;
proc sort data=Add333; by x2 x1;
proc sort data=Add474; by x2 x1;
data Table1903;
   merge Check Add333(rename=(w=w333))
               Add474(rename=(w=w474));
   by x2 x1;
   r333 = round(5*w333);
   r474 = round(9*w474);
   _rCheck = rCheck; if (_rCheck = .) then _rCheck = 0;
   _r333   = r333  ; if (_r333   = .) then _r333   = 0; t333 = _rCheck + _r333;
   _r474   = r474  ; if (_r474   = .) then _r474   = 0; t474 = _rCheck + _r474;
   drop _rCheck _r333 _r474;
   label rCheck = "Check first-order"
         r333   = "[5*w]"
         r474   = "[9*w]";
run;

title2 "Table 19.3 - Augmenting design for checking first-order model";
proc print data=Table1903 label noobs;
   var x1 x2 rCheck w333 r333 t333 w474 r474 t474;
   format w333 w474 6.4;
run;
title2;

title1;