/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 21.6                                              */
/*   TITLE: Example 21.6: Constant or quadratic model?                */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 21.6: Constant or quadratic model?";

/*
/  Candidate set
/---------------------------------------------------------------------*/
data Can;
   do x = -1 to 1 by 0.05;
      output;
      end;
run;

proc iml;
   use Can;
   read all var {x} into xCan;

   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;
      x = shape(x,nrow(x)*ncol(x),1);
      w = shape(w,nrow(w)*ncol(w),1);
      w = w/sum(w);
   finish;

   start MakeZ(x);
      F = j(nrow(x),1) || x || x#x;
      return(F);
   finish;

   %let t10 =     1;
   %let t11 =   0.5;
   %let t12 = 0.866;

   t1 = {&t10 &t11 &t12}`;

/*
/  T-optimality criterion
/---------------------------------------------------------------------*/
   start TCritWX(wx) global(t1);
      n = nrow(wx)*ncol(wx)/2;
      w = shape(wx,2,n)[1,]`;
      w = (w##2)/sum(w##2);
      x = shape(wx,2,n)[2,]`;
	  F1 = MakeZ(x); F2 = F1[,1];
      eta1 = F1*t1;
      eta2 = F2*inv(F2`*diag(w)*F2)*F2`*diag(w)*F1*t1;
      return(sum(w#((eta1 - eta2)##2)));
   finish;

   start TCritW(_w) global(x);
      w = shape(_w,nrow(_w)*ncol(_w),1);
      return(TCritWX(w` || x`));
   finish;

/*
/  D-optimality criterion for F1
/---------------------------------------------------------------------*/
   start DCritWX(wx) global(t1,eps);
      n = nrow(wx)*ncol(wx)/2;
      w = shape(wx,2,n)[1,]`;
      w = (w##2)/sum(w##2);
      x = shape(wx,2,n)[2,]`;
	  F1 = MakeZ(x);
      Dw = diag(w);
      m11 = F1`*Dw*F1;
      return(det(M11)**(1/ncol(F1)));
   finish;

   start DCritW(_w) global(x);
      w = shape(_w,nrow(_w)*ncol(_w),1);
      return(DCritWX(w` || x`));
   finish;

/*
/  DT-optimality criterion
/---------------------------------------------------------------------*/
   start DTCritW(w) global(k);
      D = DCritW(w);
      T = TCritW(w);
      if (D <= 0) then lD = -9999; else lD = log(D);
      if (T <= 0) then lT = -9999; else lT = log(T);
      return(k*lD + (1-k)*lT);
   finish;

   start DTCritWX(wx) global(k);
      D = DCritWX(wx);
      T = TCritWX(wx);
      if (D <= 0) then lD = -9999; else lD = log(D);
      if (T <= 0) then lT = -9999; else lT = log(T);
      return(k*lD + (1-k)*lT);
   finish;

   start DTOptimize;
      x = XCan;
      winit = j(nrow(X),1);
      F = MakeZ(x);
      call nlpqn(rc,w,"DTCritW",winit,1);
      wOp = w;
      wOp = ((wOp##2)/sum(wOp##2))`;
      XOp = X  [loc(wOp > 1e-6),];
      wOp = wOp[loc(wOp > 1e-6),]; wOp = wOp /sum(wOp);

      n = nrow(wOp)*ncol(wOp);

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

         wxinit = sqrt(wOp)` || xOp`;
         con =    (shape(.,1,n) || shape(-1,1,n))
               // (shape(.,1,n) || shape( 1,1,n));
         call nlpqn(rc,wxOp,"DTCritWX",wxinit,1,con);
         wOp = shape(wxOp,2,n)[1,];
         wOp = (wOp##2)/sum(wOp##2);
         xOp = shape(wxOp,2,n)[2,];
         call Consolidate(xOp,wOp,1e-2,1e-4);
         n = nrow(wOp)*ncol(wOp);

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

/*
/  Component optimum designs.
/---------------------------------------------------------------------*/
k = 0; call DTOptimize; xTOp = xOp; wTOp = wOp;   /* T-optimum design */
k = 1; call DTOptimize; xDOp = xOp; wDOp = wOp;   /* D-optimum design */

/*
/  DT-optimum designs for various values of k.
/---------------------------------------------------------------------*/
   options nonotes;
   free All;
   do k = 0 to 1 by 0.1;
      call DTOptimize;
      D = DCritWX(sqrt(wOp )` || xOp` ) / DCritWX(sqrt(wDOp)` || xDOp`);
      T = TCritWX(sqrt(wOp )` || xOp` ) / TCritWX(sqrt(wTOp)` || xTOp`);
      DT = DTCritWX(sqrt(wOp )` || xOp` );
      This = shape(k,nrow(wOp),1) || wOp || xOp || shape(D,nrow(wOp),1) || shape(T,nrow(wOp),1) || shape(DT,nrow(wOp),1);
      All = All // This;
      end;
   create All var {k w x D T DT};
   append from All;
   close All;

/*
/  Add left end-point with zero weight to first design.
/---------------------------------------------------------------------*/
data One; set All(obs=1);
   k = 0;
   w = 0;
   x = -1;
   output;
data All; set One All;
run;

proc transpose data=All out=Points (rename=(col1=x1 col2=x2 col3=x3) drop=_NAME_);
   by k;
   var x;
proc transpose data=All out=Weights(rename=(col1=w1 col2=w2 col3=w3) drop=_NAME_);
   by k;
   var w;
proc transpose data=All out=DEff(rename=(col1=DEff) drop=_NAME_ col2 col3);
   by k;
   var D;
proc transpose data=All out=TEff(rename=(col1=TEff) drop=_NAME_ col2 col3);
   by k;
   var T;
data fig21DT; merge Points Weights DEff TEff;
run;

title2 "Fig. 21.5 -  Structure of DT-optimum designs as k varies";
%let pagesize = %sysfunc(getoption(pagesize));
options ps=40;
proc plot data=fig21DT;
   plot x1*k='1' x2*k='2' x3*k='3' / overlay;
   plot w1*k='1' w2*k='2' w3*k='3' / overlay;
   plot DEff*k='D' TEff*k='T'      / overlay;
run;
options ps=&pagesize;
title2;

title1;