/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 12.2                                              */
/*   TITLE: Equation 12.2 - Exact design weights over [-1,0,1]**2     */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Equation 12.2 - Exact design weights over [-1,0,1]**2";



%let nInc = 11;
proc factex;
   factors x1 x2 / nlev = &nInc;
   output out=Grid;
data Grid; set Grid;
   x1 = 2*x1/(&nInc-1) - 1;
   x2 = 2*x2/(&nInc-1) - 1;
run;


proc iml;

   /*
   /  Define the basic criterion as a function of both the weight and
   /  the support of the design.  For this problem, the criterion is
   /  the variance of the second parameter, the one corresponding to
   /  the x*x term.  This is the [2,2] element of the variance matrix
   /  for the parameters.
   /------------------------------------------------------------------*/
   start Crit(w,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`*diag(w)*Z/nrow(x);
      return(log(det(Mi))/ncol(Mi));
   finish;

   start CritW(w) global(x);            /* Crit for sqrt(w), x fixed */
      return(Crit((w##2)/sum(w##2),x));
   finish;

   /*
   /  This function simplifies a design, often revealing an equivalent
   /  design with fewer support and thus one whose structure is easier
   /  to see.
   /------------------------------------------------------------------*/
   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 Consolidate(x,w,ex,ew);

      x = x[loc(w > ew),];             /* Drop points with small wgts */
      w = w[loc(w > ew)];

      Group = shape(0,1,ncol(x));      /* Replace nearby points with  */
      do i = 1 to ncol(x);             /* their avg, with wgts equal  */
         if (^Group[i]) then do;       /* to sum(wgt) over averaged   */
            Group[i] = max(Group) + 1; /* points.                     */
            do j = i+1 to ncol(x);
               if (max(abs(x[i,] - x[j,])) < ex) then
                  Group[j] = Group[i];
               end;
            end;
         end;
      xNew = shape(.,ncol(x),max(Group));
      wNew = shape(.,1      ,max(Group));
      do i = 1 to max(Group);
         xNew[i,] = (x[loc(Group = i),])[,:];
         wNew[i ] = (w[loc(Group = i) ])[ +];
         end;

      wx = w || x;
      call sort(wx,1:ncol(x)+1);       /* Make sure the points are    */
      w = wx[,1];                      /* sorted by w & x.            */
      x = wx[,2:ncol(x)+1];
   finish;

   /*
   /  A little function to turn a number into a string for printing.
   /------------------------------------------------------------------*/
   start tlc(n); return(trim(left(char(n,12)))); finish;


   use Grid; read all var {"x1" "x2"} into xGrid;

   x = xGrid;
   winit = j(nrow(x),1) / nrow(x);
   winit = sqrt(winit);
   call nlpqn(rc,wOp,"CritW",winit,1);
   wOp = (wOp##2)`/sum(wOp##2);
   xOp = x;
   print "Optimal weights over whole grid",,wOp [format=7.4] xOp;

   call Consolidate(xOp,wOp,0,1e-4);
   print "Consolidated design",,wOp [format=7.4] xOp;
   x = xOp;
   winit = sqrt(wOp);
   call nlpqn(rc,wOp,"CritW",winit,1);
   wOp = (wOp##2)`/sum(wOp##2);
   xOp = x;
   print "Improved consolidated design",,wOp [format=7.4] xOp;

title1;