/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 9.1                                               */
/*   TITLE: Example 9.1 - Quadratic regression through the origin     */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: This is a complicated program because of its generality:  */
/*          the same code can be used with minimal changes for        */
/*          finding any optimal design for a single factor, and the   */
/*          approach is also useful for designs with more factors.    */
/*          IML functions are defined to compute the criterion for    */
/*          both w and x, for w with x fixed, and for x with w fixed, */
/*          and these functions are iteratively optimized to find     */
/*          both the optimum w and the optimum x.                     */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Computations for Figure 6.9 - Variance dispersion graph";


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 = x || x#x;
      Mi = Z`*diag(w)*Z;
      return(inv(Mi)[2,2]);
   finish;

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

   start CritX(x) global(w);            /* Crit for      x , w fixed */
      return(Crit(w,x));
   finish;

   start CritWX(wx);                    /* Crit for sqrt(w) & x      */
      w = shape(wx,2,ncol(wx)/2)[1,]`;
      w = (w##2)/sum(w##2);
      x = shape(wx,2,ncol(wx)/2)[2,]`;
      return(Crit(w,x));
   finish;

   /*
   /  This function simplifies a design, often revealing an equivalent
   /  design with fewer support and thus one whose structure is easier
   /  to see.
   /------------------------------------------------------------------*/
   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 (abs(x[i] - x[j]) < ex) then
                  Group[j] = Group[i];
               end;
            end;
         end;
      xNew = shape(.,1,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;

      r = rank(xNew);                  /* Make sure the points are    */
      x = xNew;                        /* sorted by x.                */
      w = wNew;
      x[r] = xNew;
      w[r] = wNew;
   finish;

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

/*
/  Initial optimization: Fix the support at a grid values between 0
/  and 1 and optimize w.
/---------------------------------------------------------------------*/
%let xLo = 0;
%let xHi = 1;
%let Inc = 40;

   x    = &xLo + (&xHi - &xLo)*(0:&Inc)`/&Inc;
   winit = j(nrow(x),1) / nrow(x);
   winit = sqrt(winit);
   call nlpqn(rc,wOp,"CritW",winit,0);
   wOp = (wOp##2)/sum(wOp##2);
   xOp = x`;

   print "Initial optimum weights over support grid",,
         (xOp`)            [colname={"Support"}  format=7.4]
         (wOp`)            [colname={"Weights"}  format=7.4]
         (Crit(wOp`,xOp`)) [colname={"Var[2,2]"} format=6.4];

/*
/  Most of the weights in the inial optimization are near zero.
/  Simplifiy the design and redisplay.
/---------------------------------------------------------------------*/
   call Consolidate(xOp,wOp,0,1e-4);
   print "Simplified inital optimum design",,
         (xOp`)            [colname={"Support"}  format=7.4]
         (wOp`)            [colname={"Weights"}  format=7.4]
         (Crit(wOp`,xOp`)) [colname={"Var[2,2]"} format=6.4];

/*
/  Now, starting from this design, iterate through jointly optimizing
/  w & x and simplifying the result.  In this case, the process
/  converges after just one joint-optimization, but in other cases
/  it can take longer for the structure of the optimum design to
/  become manifest
/---------------------------------------------------------------------*/
   i = 0;
   Improvement = 1;
   do while ((i < 100) & (Improvement > 1e-4));
      i = i + 1;

      last = (xOp` || wOp`);

      wxinit = sqrt(wOp) || xOp;
      con =    (shape(.,1,ncol(xOp)) || shape(&xLo,1,ncol(xOp)))
            // (shape(.,1,ncol(xOp)) || shape(&xHi,1,ncol(xOp)));
      call nlpqn(rc,wxOp,"CritWX",wxinit,0,con);
      wOp = shape(wxOp,2,ncol(wxOp)/2)[1,];
      xOp = shape(wxOp,2,ncol(wxOp)/2)[2,];
      wOp = (wOp##2)/sum(wOp##2);
      call Consolidate(xOp,wOp,1e-2,1e-5);

      if (ncol(xOp) ^= nrow(last)) then
         Improvement = 999;
      else
         Improvement = max(abs(last - (xOp` || wOp`)));

      print ("Improving both weight and support: iteration " + tlc(i)),,
         (xOp`)            [colname={"Support"}  format=7.4]
         (wOp`)            [colname={"Weights"}  format=7.4]
         (Crit(wOp`,xOp`)) [colname={"Var[2,2]"} format=6.4]
         Improvement;
      end;

/*
/  Display the final optimum.
/---------------------------------------------------------------------*/
   print "Final optimum design",,
         (xOp`)            [colname={"Support"}  format=7.4]
         (wOp`)            [colname={"Weights"}  format=7.4]
         (Crit(wOp`,xOp`)) [colname={"Var[2,2]"} format=6.4];
quit;
title1;