/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 13.2                                              */
/*   TITLE: SAS Task 13.2 - Exact optimal quadratic designs on square */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Same as Program 12.1 - Example 12.1                       */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "SAS Task 13.2 - Exact optimal quadratic designs on square";


/*
/  Designs of Box and Draper (1971)
/---------------------------------------------------------------------*/
proc factex;
   factors x1 x2;
   output out=Square;
run;

data BD6;
   x1 =  -1; x2 =  -1; output;
   x1 =   1; x2 =  -1; output;
   x1 =  -1; x2 =   1; output;
   a = (4 - sqrt(13))/3;
   x1 =  -a; x2 =  -a; output;
   x1 =   1; x2 = 3*a; output;
   x1 = 3*a; x2 =   1; output;
run;

data BD7;
   x1 = -0.092; x2 =  0.092; output;
   x1 =  1    ; x2 = -0.067; output;
   x1 =  0.067; x2 = -1    ; output;
data BD7; set BD7 Square;
run;

data BD8;
   x1 =  1    ; x2 =  0    ; output;
   x1 =  0.082; x2 =  1    ; output;
   x1 =  0.082; x2 = -1    ; output;
   x1 = -0.215; x2 =  0    ; output;
data BD8; set BD8 Square;
run;

proc factex;
   factors x1 x2 / nlev=3;
   output out=BD9;
run;
   


/*
/  Use constrained nonlinear optimization to find optimal designs
/  directly.  Begin by defining a function to compute the log
/  determininant for a quadratic model.
/---------------------------------------------------------------------*/
options nonotes;
proc iml;

   start llog(x); if (x > 0) then return(log(x)); else return(-9999); finish;
   start Crit(x);
      Z = j(nrow(x),1) || x;
      do i = 1 to ncol(x);
         do j = i to ncol(x);
            Z = Z || x[,i]#x[,j];
            end;
         end;
      Mi = Z`*Z/nrow(x);
      return(llog(det(Mi))/ncol(Mi));
   finish;
   start CritX(_x) global(m,N);
      return(Crit(shape(_x,N,m)));
   finish;

   start tlc(n); return(trim(left(char(n,12)))); finish;

/*
/  A function to find the best exact design for a given number of
/  factors and points.  The approach uses 10 random restarts; this
/  might have to be modified for more than 2 factors.
/---------------------------------------------------------------------*/
   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 OpDesign(_m,_N) global(N,m);
      N = _N;
      m = _m;
      BestD  = -9999;
      free BestX;
      do i = 1 to 10;
         xinit = 2*ranuni(1:(N*m))-1;
         con =    shape(-1,1,N*m)
               // shape( 1,1,N*m);
         call nlpqn(rc,xOp,"CritX",xinit,1,con);
         if (CritX(xOp) > BestD) then do;
            BestD = CritX(xOp);
            BestX = xOp;
            end;
         end;
      x = shape(BestX,N,m);
      call sort(x,1:ncol(x));
      return(x);
   finish;

/*
/  Read in Box and Draper's designs, compute corresponding designs
/  directly, and compare.
/---------------------------------------------------------------------*/

   use BD6; read all var {x1 x2} into BD6; close;
   use BD7; read all var {x1 x2} into BD7; close;
   use BD8; read all var {x1 x2} into BD8; close;
   use BD9; read all var {x1 x2} into BD9; close;

   Op6 = OpDesign(2,6);
   Op7 = OpDesign(2,7);
   Op8 = OpDesign(2,8);
   Op9 = OpDesign(2,9);

/*
/  The results confirm that Box and Draper's designs are D-optimal.
/---------------------------------------------------------------------*/
   print "Optimal designs and relative efficiency of BD's",,
         (1*Op6) [format=7.4] (exp(Crit(BD6))/exp(Crit(Op6))) [format=7.4],,
         (1*Op7) [format=7.4] (exp(Crit(BD7))/exp(Crit(Op7))) [format=7.4],,
         (1*Op8) [format=7.4] (exp(Crit(BD8))/exp(Crit(Op8))) [format=7.4],,
         (1*Op9) [format=7.4] (exp(Crit(BD9))/exp(Crit(Op9))) [format=7.4];

title1;