/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 11.3                                              */
/*   TITLE: Table 11.2 - Efficiencies for Example 1.1                 */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Table 11.2 - Efficiencies for Example 1.1";


/*
/  Data for Table 1.1 - Desorption of carbon monoxide
/---------------------------------------------------------------------*/
data CO;
   input KCRatio n @@;
   do i = 1 to n;
      input Desorbed @@;
      output;
      end;
datalines;
0.05 2 0.05 0.10
0.25 2 0.25 0.35
0.50 3 0.75 0.85 0.95
1.25 5 1.42 1.75 1.82 1.95 2.45
2.10 6 3.05 3.19 3.25 3.43 3.50 3.93
2.50 4 3.75 3.93 3.99 4.07
;


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

options nonotes;
proc iml;
   use CO; read all var {KCRatio} into xCO;
   xCO = xCO/max(xCO);

   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 Evaluate(CritFunction,Label);
   start Crit(w,x); return(&CritFunction(w,x)); finish;

%let xLo = 0;
%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`;
   DOpt = exp(Crit(wOp,xOp));
   DCO  = exp(Crit(j(nrow(xCO),1)/nrow(xCO),xCO));
   print &Label,,
         (1*xOp)    [colname={"Support"}          format=7.4]
         (1*wOp)    [colname={"Weights"}          format=7.4]
         (DCO/DOpt) [colname={"Efficiency of CO"} format=6.4];
%mend;

   start Crit1(w,x);
      Z = j(nrow(x),1) || x;
      Mi = Z`*diag(w)*Z;
      return(log(det(Mi))/ncol(Mi));
   finish;

   start Crit2(w,x);
      Z = j(nrow(x),1) || x || x#x;
      Mi = Z`*diag(w)*Z;
      return(log(det(Mi))/ncol(Mi));
   finish;

   start Crit3(w,x);
      Z = j(nrow(x),1) || x || x#x;
      Mi = Z`*diag(w)*Z;
      return(-log(inv(Mi)[3,3]));
   finish;

   start Crit4(w,x);
      Z =  x;
      Mi = Z`*diag(w)*Z;
      return(log(det(Mi))/ncol(Mi));
   finish;

   start Crit5(w,x);
      Z = x || x#x;
      Mi = Z`*diag(w)*Z;
      return(log(det(Mi))/ncol(Mi));
   finish;

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

%Evaluate(Crit1,"Final optimum linear design");
%Evaluate(Crit2,"Final optimum quadratic design");
%Evaluate(Crit3,"Final optimum design for b2 in quadratic model");
%Evaluate(Crit4,"Final optimum linear design without intercept");
%Evaluate(Crit5,"Final optimum quadratic design without intercept");
%Evaluate(Crit6,"Final optimum design for b2 in quadratic model without intercept");