/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 21.1                                              */
/*   TITLE: Section 21.4 and SAS Task 21.1 - Compound optimum designs */
/*          for polynomials in one factor                             */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Compound optimum designs for polynomials in one factor";

/*
/  All very much the same as other programs for computing continuous
/  optimum designs using nonlinear optimization, except that now the
/  function which makes the design matrix depends on the order of the
/  polynomial, and the critical value function DCrit() depends on the
/  vector of compounding parameters kappa.
/---------------------------------------------------------------------*/
proc iml;

   start MakeZ(x,k);
      free F;
      term = j(nrow(x),1);
      do d = 0 to k;
         F = F || term;
         term = term#x;
         end;
      return(F);
   finish;

   start DCrit(w,x) global(kappa);
      dd = 0;
      do d = 2 to 6;
         Fd = MakeZ(x,d);
         dd = dd + kappa[d]*llog(det(Fd`*diag(w)*Fd))/ncol(Fd);
         end;
      return(dd);
   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 llog(x); if (x > 0) then return(log(x)); else return(-9999); finish;

   start DCritW(w) global(x);
      return(DCrit((w##2)/sum(w##2),x));
   finish;

   start DCritX(x) global(w);
      return(DCrit(w,x`));
   finish;
   
   start DCritWX(wx);
      w = shape(wx,2,ncol(wx)/2)[1,];
      w = (w##2)/sum(w##2);
      x = shape(wx,2,ncol(wx)/2)[2,];
      return(DCrit(w`,x`));
   finish;

   start Optimize;
      nint = 51;
      x    = -1 + (1 - -1)*((1:nint)`-1)/(nint-1);
      winit = j(nrow(x),1) / nrow(x);
      call nlpqn(rc,wOp,"DCritW",winit,1);
      wOp = ((wOp##2)/sum(wOp##2))`;
      xOp = x  [loc(wOp > 1e-5)]`;
      wOp = wOp[loc(wOp > 1e-5)]`; wOp = wOp /sum(wOp);

      i = 0;
      diff = 1;
      do while ((i < 100) & (diff > 1e-6));
         i = i + 1;
         last = (xOp` || wOp`);

         wxinit = sqrt(wOp) || xOp;
         con =    (shape(.,1,ncol(xOp)) || shape(-1,1,ncol(xOp)))
               // (shape(.,1,ncol(xOp)) || shape( 1,1,ncol(xOp)));
         call nlpqn(rc,wxOp,"DCritWX",wxinit,1,con);
         wOp = shape(wxOp,2,ncol(wxOp)/2)[1,];
         wOp = (wOp##2)/sum(wOp##2);
         xOp = shape(wxOp,2,ncol(wxOp)/2)[2,];
         call Consolidate(xOp,wOp,1e-2,1e-3);

         if (ncol(xOp) ^= nrow(last)) then
            diff = 1;
         else
            diff = max(abs(last - (xOp` || wOp`)));
         end;
   xOp = xOp`;
   wOp = wOp`;
   finish;

/*
/  Compute each individual optimum design for polynomials of
/  orders 2 to 6.
/---------------------------------------------------------------------*/
   kappa = ((1:6)  = 2); run Optimize; wOp2 = wOp; xOp2 = xOp;
   kappa = ((1:6)  = 3); run Optimize; wOp3 = wOp; xOp3 = xOp;
   kappa = ((1:6)  = 4); run Optimize; wOp4 = wOp; xOp4 = xOp;
   kappa = ((1:6)  = 5); run Optimize; wOp5 = wOp; xOp5 = xOp;
   kappa = ((1:6)  = 6); run Optimize; wOp6 = wOp; xOp6 = xOp;

/*
/  Compute the compound optimum design.
/---------------------------------------------------------------------*/
   kappa = ((1:6) <= 6); run Optimize; wOpC = wOp; xOpC = xOp;

   print " DOpt 6th Order" "  " "Compound DOpt",
         "   Wgt   Supp" "   " "   Wgt   Supp",,
         (wOp6 || xOp6) [format=6.4] "   "
         (wOpC || xOpC) [format=6.4];

/*
/  Compute efficiencies of D-optimum 6th order and compound D-optimum
/  designs.
/---------------------------------------------------------------------*/
   kappa = ((1:6)  = 2); DEff26 = exp(DCrit(wOp6,xOp6))/exp(DCrit(wOp2,xOp2));
                         DEff2C = exp(DCrit(wOpC,xOpC))/exp(DCrit(wOp2,xOp2));

   kappa = ((1:6)  = 3); DEff36 = exp(DCrit(wOp6,xOp6))/exp(DCrit(wOp3,xOp3));
                         DEff3C = exp(DCrit(wOpC,xOpC))/exp(DCrit(wOp3,xOp3));

   kappa = ((1:6)  = 4); DEff46 = exp(DCrit(wOp6,xOp6))/exp(DCrit(wOp4,xOp4));
                         DEff4C = exp(DCrit(wOpC,xOpC))/exp(DCrit(wOp4,xOp4));

   kappa = ((1:6)  = 5); DEff56 = exp(DCrit(wOp6,xOp6))/exp(DCrit(wOp5,xOp5));
                         DEff5C = exp(DCrit(wOpC,xOpC))/exp(DCrit(wOp5,xOp5));

   kappa = ((1:6)  = 6); DEff66 = exp(DCrit(wOp6,xOp6))/exp(DCrit(wOp6,xOp6));
                         DEff6C = exp(DCrit(wOpC,xOpC))/exp(DCrit(wOp6,xOp6));

   print "Table 21.1 - D-efficiencies of D-optimum and compound D-optimum designs",,
         (  (DEff26 || DEff2C)
          //(DEff36 || DEff3C)
          //(DEff46 || DEff4C)
          //(DEff56 || DEff5C)
          //(DEff66 || DEff6C)) [format=percent9.2];

title1;