/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 11.4                                              */
/*   TITLE: Table 11.3 - Polynomial regression in one variable        */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Table 11.3 - Polynomial regression in one variable";


/*
/  The continuous optimization proceeds very much as in Program 9.1,
/  except for the different models.
/---------------------------------------------------------------------*/

options nonotes;
proc iml;

   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;



%macro PolyOpt(k);
   k = &k;

%let xLo = -1;
%let xHi =  1;
%let Inc = 40;

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

   call Consolidate(xOp,wOp,0,1e-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,1,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`)));
      end;
   
   wOp = wOp`;
   xOp = xOp`;
%mend;


   start Crit(w,x) global(k);
      Z = j(nrow(x),1);
      xi =  j(nrow(x),1);
      do i = 1 to k;
         xi = xi#x;
         Z = Z || xi;
         end;
      Mi = Z`*diag(w)*Z;
      return(log(det(Mi))/ncol(Mi));
   finish;


%PolyOpt(1); Support = xOp`;
             Weights = wOp`;
%PolyOpt(2); Support = (Support || shape(.,nrow(Support),1)) // xOp`;
             Weights = (Weights || shape(.,nrow(Weights),1)) // wOp`;
%PolyOpt(3); Support = (Support || shape(.,nrow(Support),1)) // xOp`;
             Weights = (Weights || shape(.,nrow(Weights),1)) // wOp`;
%PolyOpt(4); Support = (Support || shape(.,nrow(Support),1)) // xOp`;
             Weights = (Weights || shape(.,nrow(Weights),1)) // wOp`;
%PolyOpt(5); Support = (Support || shape(.,nrow(Support),1)) // xOp`;
             Weights = (Weights || shape(.,nrow(Weights),1)) // wOp`;
%PolyOpt(6); Support = (Support || shape(.,nrow(Support),1)) // xOp`;
             Weights = (Weights || shape(.,nrow(Weights),1)) // wOp`;
%PolyOpt(7); Support = (Support || shape(.,nrow(Support),1)) // xOp`;
             Weights = (Weights || shape(.,nrow(Weights),1)) // wOp`;

Order = (1:7)`;
print Order Support [format=7.4],,
      Order Weights [format=7.4];