/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 10.1                                              */
/*   TITLE: Figure 10.1 - Quadratic regression through the origin     */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: This program is largely the same as Program 9.1, which    */
/*          see for comments on the method.  The computation at the   */
/*          end of the program of the variance function for the Ds    */
/*          criterion is the only difference.                         */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Figure 10.1 - Quadratic regression through the origin";

proc iml;

   start Crit(w,x);
      Z = x || x#x;
      Mi = Z`*diag(w)*Z;
      return(inv(Mi)[2,2]);
   finish;

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

   start CritX(x) global(w);
      return(Crit(w,x));
   finish;

   start CritWX(wx);
      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;

   start Consolidate(x,w,ex,ew);

      x = x[loc(w > ew)]`;
      w = w[loc(w > ew)]`;

      Group = shape(0,1,ncol(x));
      do i = 1 to ncol(x);
         if (^Group[i]) then do;
            Group[i] = max(Group) + 1;
            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);
      x = xNew;
      w = wNew;
      x[r] = xNew;
      w[r] = wNew;
   finish;

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


%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];

   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];

   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;
   
   wOp = wOp`;
   xOp = xOp`;
   print "Final optimum design",,
         (1*xOp)         [colname={"Support"}  format=7.4]
         (1*wOp)         [colname={"Weights"}  format=7.4]
         (Crit(wOp,xOp)) [colname={"Var[2,2]"} format=6.4];

   f1 = xOp#xOp;
   f2 = xOp;
   F = f1 || f2;
   M = F`*diag(wOp)*F;

   start d(x) global(M);
      f1 = x#x;
      f2 = x;
      f = f1 || f2;
      return(F*inv(M)*F` - f2`*inv(M[2,2])*f2);
   finish;

   free dsx;
   do x = 0 to 1 by 0.1;
      dsx = dsx // (x  || d(x));
      end;
   create dsx var {"x" "d"}; append from dsx;
quit;

title2 "Variance function ds(x,xi)";
proc print data=dsx;
run;
title2;
title1;