/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 14.1                                              */
/*   TITLE: Example 14.1 - Second-order response surface              */
/*          with one qualitative factor                               */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 14.1 - Second-order response surface with one qualitative factor";

/*
/  Create 3**3 candidate set and use OPTEX to generate the optimal
/  blocked design.
/---------------------------------------------------------------------*/
title2 "Optimal 9-run 3-block quadratic design over 3**2 grid";
data GCandidates1;
      do Group = 1 to 3;
         do x1 = -1 to 1;
            do x2 = -1 to 1;
               output;
               end;
            end;
         end;
proc optex data=GCandidates1 coding=orthcan;
   class Group;
   model Group x1 x2 x1*x2 x1*x1 x2*x2;
   generate n=9 method=m_fedorov niter=1000 keep=10;
   output out=Design9;
proc sort data=Design9;
   by x1 x2 Group;
proc print data=Design9;
run;
title2;

/*
/  Use optimal blocking in OPTEX with INITDESIGN= and NOEXCHANGE in
/  order to optimally block a fixed design for x1 and x2.
/---------------------------------------------------------------------*/
title2 "Optimally blocked 3**2 design";
data Factorial;
   do x1 = -1 to 1;
      do x2 = -1 to 1;
         output;
         end;
      end;
proc optex data=Factorial coding=orthcan;
   model x1 x2 x1*x2 x1*x1 x2*x2;
   generate initdesign=Factorial method=sequential;
   block structure=(3)3 noexchange niter=100 keep=10;
   output out=GFactorial(rename=(Block=Group));
proc sort data=GFactorial;
   by x1 x2 Group;
proc print data=GFactorial;
run;
title2;

/*
/  Create the grid over which d_ave and d_max will be computed.
/---------------------------------------------------------------------*/
data Grid;
   do Group = 1 to 3;
      do x1 = -1 to 1 by 0.1;
         do x2 = -1 to 1 by 0.1;
            output;
            end;
         end;
      end;
run;


proc iml;
   /*
   /  Read in points of the optimally grouped factorial design, ...
   /------------------------------------------------------------------*/
   use GFactorial; read all var ("x1":"x2"); read all var ("Group");

   /*
   /  ... use them to create the design matrix for this design, ...
   /------------------------------------------------------------------*/
   XF = design(Group) || x1 || x2 || x1#x2 || x1#x1 || x2#x2;

   /*
   /  ... and then the information.
   /------------------------------------------------------------------*/
   MF = XF`*XF / nrow(XF);

   /*
   /  Do the same for the optimal design.
   /------------------------------------------------------------------*/
   use Design9; read all var ("x1":"x2"); read all var ("Group");
   XD = design(Group) || x1 || x2 || x1#x2 || x1#x1 || x2#x2;
   MD = XD`*XD / nrow(XD);

   p = ncol(XD);

   /*
   /  Print the determinants and the ratio of their p-throots; note
   /  that this fairly close to the ratio of the numbers in AD.
   /------------------------------------------------------------------*/
   detGrp = det(MF);
   detOpt = det(MD);
   RelEff = (detGrp**(1/p))/(detOpt**(1/p));
   print "D for optimal and optimally grouped designs, and relative efficiency",,
      detOpt detGrp ((detGrp**(1/p))/(detOpt**(1/p))) [format=percent9.2];


   /*
   /  Read in the points over which d_ave and d_max will be computed
   /  and create their design matrix.
   /------------------------------------------------------------------*/
   use Grid;
   read all var ("x1":"x2"); read all var ("Group");
   XE = design(Group) || x1 || x2 || x1#x2 || x1#x1 || x2#x2;

   /*
   /  Compute the variances of prediction.
   /------------------------------------------------------------------*/
   free dXD dXF;
   do i = 1 to nrow(XE);
      dXD = dXD // XE[i,]*inv(MD)*XE[i,]`;
      dXF = dXF // XE[i,]*inv(MF)*XE[i,]`;
      end;

   dXD = vecdiag(XE*inv(MD)*XE`);        /* Optimal design            */
   dXF = vecdiag(XE*inv(MF)*XE`);        /* Grouped factorial design  */

   /*
   /  And this is matrix-programming shorthand for the average and
   /  the maximum, respectively, of these two vectors.
   /------------------------------------------------------------------*/
   dAve = dXD[: ] // dXF[: ];
   dMax = dXD[<>] // dXF[<>];
   print ("Optimal" // "Grouped") dAve [format=5.2] dMax [format=5.2];


/*
/  Create 3*[-1,1]**2 candidate set and use OPTEX to generate the
/  optimal blocked design.
/---------------------------------------------------------------------*/
title2 "Optimal 13-run 2-block quadratic design over 3**2 grid";
data GCandidates2;
      do Group = 1 to 2;
         do x1 = -1 to 1;
            do x2 = -1 to 1;
               output;
               end;
            end;
         end;
proc optex data=GCandidates2 coding=orthcan;
   class Group;
   model Group x1 x2 x1*x2 x1*x1 x2*x2;
   generate n=13 method=m_fedorov niter=1000 keep=10;
   output out=Design13;
proc sort data=Design13;
   by Group x1 x2;
proc print data=Design13;
   by Group;
run;
title2;


title2 "Improved 9-run 3-block quadratic design over [-1,1]**2";
proc iml;
   start MakeZ(B,x);
      Z = design(B);
      do i = 1 to ncol(x);
         Z = Z || x[,i];
         end;
      do i = 1 to ncol(x)-1;
         do j = i+1 to ncol(x);
            Z = Z || x[,i]#x[,j];
            end;
         end;
      do i = 1 to ncol(x);
         Z = Z || x[,i]#x[,i];
         end;
      return(Z);
   finish;

   start DCrit(xx) global(B);
      x = shape(xx,nrow(xx)*ncol(xx)/2,2);
      Z = MakeZ(B,x);
      d = det(Z`*Z/nrow(Z));
      if (d > 0) then return(log(d)/ncol(Z));
      else            return(-9999);
   finish;

   use Design9;
   read all var ("x1":"x2") into xop;
   read all var ("Group")   into B;

   /*
   /  Optimally adjust points of discrete optimal design.
   /------------------------------------------------------------------*/
   xinit = shape(xop,2*nrow(xop),1);
   dOpDesign = DCrit(xinit);
   con =    shape(-1,1,nrow(xinit))
         // shape( 1,1,nrow(xinit));
   call nlpqn(rc,xx,"DCrit",xinit,1,con);
   dxx = DCrit(xx);
   xim = shape(xx,nrow(xx)*ncol(xx)/2,2);

   /*
   /  Evaluate dCrit for design measure originally claimed to be
   /  optimum.
   /------------------------------------------------------------------*/
   use GCandidates1;
   read all var ("x1":"x2") into xCan;
   read all var ("Group")   into BCan;

   Corners     = loc(    (abs(xCan[,1])=1) & (abs(xCan[,2])=1) );
   EdgeCenters = loc(  ^((abs(xCan[,1])=1) & (abs(xCan[,2])=1))
                     &  ((abs(xCan[,1])=1) | (abs(xCan[,2])=1)));
   Center      = loc(    (abs(xCan[,1])=0) & (abs(xCan[,2])=0) );
   w = shape(.,nrow(xCan),1);
   w[Corners    ] = 0.1458;
   w[EdgeCenters] = 0.0802;
   w[Center     ] = 0.0960;
   w = w / sum(w);
   Z = MakeZ(BCan,xCan);
   dCan1 = log(det(Z`*diag(w)*Z))/ncol(Z);

   /*
   /  Optimize dCrit over the design measure.
   /------------------------------------------------------------------*/
   start DCritW(w) global(Z);
      Dw = diag((w##2)/sum(w##2));
      return(log(det(Z`*Dw*Z))/ncol(Z));
   finish;

   winit = shape(1,nrow(Z),1); w = w / sum(w);
   call nlpqn(rc,w,"DCritW",winit,1);
   w = (w##2)/sum(w##2);

   /*
   /  Print optimal weights on three different groups of points.
   /------------------------------------------------------------------*/
   isCorner     = (    (abs(xCan[,1])=1) & (abs(xCan[,2])=1) );
   isEdgeCenter = (  ^((abs(xCan[,1])=1) & (abs(xCan[,2])=1))
                   &  ((abs(xCan[,1])=1) | (abs(xCan[,2])=1)));
   isCenter     = (    (abs(xCan[,1])=0) & (abs(xCan[,2])=0) );
   print "Optimal weights over 3 * 3**2 grid",,
         xCan (w`) isCorner isEdgeCenter isCenter;
   print "Optimal weights by point type",,
         (3*(   (w[Corners    ])[:]
             // (w[EdgeCenters])[:]
             // (w[Center     ])[:])) [format=7.4];

   /*
   /  Evaluate dCrit for design measure found above to be optimum
   /  and compare with original.
   /------------------------------------------------------------------*/
   dCan2 = log(det(Z`*diag(w)*Z))/ncol(Z);

   print "D for optimal weighted design over 3 * 3**2 grid, by two methods",,
         (exp(dCan1) // exp(dCan2)) [format=best12.],
         (exp(dCan1)  / exp(dCan2)) [format=best12.];

   /*
   /  Print relative efficiencies of exact designs.
   /------------------------------------------------------------------*/
   GridDesign     = exp(dOpDesign)/exp(dCan2);
   ImprovedDesign = exp(dxx)      /exp(dCan2);
   print "D-efficiencies of design over 3 * 3**2 grid and improvement",,
         GridDesign     [format=percent8.2]
         ImprovedDesign [format=percent8.2];

   xim = B || xim;
   create imDesign var ("Block"||"x1"||"x2");
   append from xim;

proc sort  data=imDesign; by Block x1 x2;
proc print data=imDesign;
   by Block;
run;

title2;

title1;