/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 20.5                                              */
/*   TITLE: Figure 20.4 - Derivative function for T-optimum design    */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Figure 20.4 - Derivative function for T-optimum design";

%let t10 =  4.5;
%let t11 = -1.5;
%let t12 = -2.0;

/*
/  Linear terms and true model.
/---------------------------------------------------------------------*/
data Can;
   do x = -1 to 1 by 0.05;
      x10 = 1; x11 = exp( x); x12 = exp(-x);
      x20 = 1; x21 = x;       x22 = x*x;
      eta1 = &t10*x10 + &t11*x11 + &t12*x12;
      w = 0;
      output;
      end;
proc iml;
   use Can;
   read all var {x} into xCan;

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

   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 D2(w,F1,F2,t1);           /* Equations 20.21 - 20.23 in text */
      eta1 = F1*t1;
      eta2 = F2*inv(F2`*diag(w)*F2)*F2`*diag(w)*F1*t1;
      return(log(sum(w#((eta1 - eta2)##2))));
   finish;

   start D2W(sw) global(F1,F2,t1); /* Over w for given x              */
      w = shape(sw,nrow(sw)*ncol(sw),1);
      w = (w##2)/sum(w##2);
	  return(D2(w,F1,F2,t1));
   finish;

   start D2X(_x) global(w,t1);     /* Over x for given w              */
      n = nrow(_x)*ncol(_x);
      x = shape(_x,n,1);
      F1 = j(n,1) || exp(x) || exp(-x);
      F2 = j(n,1) || x      || x#x;
	  return(D2(w,F1,F2,t1));
   finish;

   start D2WX(wx) global(t1);      /* Over both w and x               */
      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 = j(n,1) || exp(x) || exp(-x);
      F2 = j(n,1) || x      || x#x;
	  return(D2(w,F1,F2,t1));
   finish;

   start D2maxWX(wx) global(t1);   /* Max equation 20.24 in text      */
      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 = j(n,1) || exp(x) || exp(-x);
      F2 = j(n,1) || x      || x#x;
      t2 = inv(F2`*diag(w)*F2)*F2`*diag(w)*F1*t1;

      psi2max = -1;
      do xx = -1 to 1 by 0.01;
         fF1 = 1 || exp(xx) || exp(-xx);
         fF2 = 1 || xx      || xx#xx;
         eta1 = fF1*t1;
         eta2 = fF2*t2;
         psi2 = (eta1 - eta2)**2;
         if (psi2 > psi2max) then psi2max = psi2;
         end;
      return(log(psi2max));
   finish;

   xOp = xCan;
   wOp = j(nrow(xOp),1);

   x = xOp;
   winit = wOp;
   F1 = j(nrow(x),1) || exp(x) || exp(-x);
   F2 = j(nrow(x),1) || x      || x#x;
   call nlpqn(rc,w,"D2W",winit,1);
   wOp = w;
   wOp = ((wOp##2)/sum(wOp##2))`;
   xOp = x;

   xOp = xOp[loc(wOp > 1e-6)];
   wOp = wOp[loc(wOp > 1e-6)]; wOp = wOp /sum(wOp);
   n = nrow(wOp)*ncol(wOp);

   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,"D2WX",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-4,1e-5);
      n = nrow(wOp)*ncol(wOp);

      if (n ^= nrow(last)) then
         diff = 1;
      else
         diff = max(abs(last - (xOp || wOp)));
      end;
   Del2 = exp(D2maxWX(sqrt(wOp)` || xOp`));
   print "Equation 20.28 - T-optimum design",,
         (xOp` // wOp`) [format=7.4] Del2;

   data = xOp || wOp;
   create DOpt var {x w}; append from data;
data DOpt; set DOpt;
   x10 = 1; x11 = exp( x); x12 = exp(-x);
   x20 = 1; x21 = x;       x22 = x*x;
   eta1 = &t10*x10 + &t11*x11 + &t12*x12;
run;

data fig2004;
   do x = -1 to 1 by 0.005;
      x10 = 1; x11 = exp( x); x12 = exp(-x);
      x20 = 1; x21 = x;       x22 = x*x;
      eta1 = &t10*x10 + &t11*x11 + &t12*x12;
      w = 0;
      output;
      end;
data DOpt; set DOpt fig2004;
proc reg data=DOpt noprint;
   model eta1 = x21 x22;
   weight w;
   output out=rDOpt r=r;
data rDOpt; set rDOpt; psi2 = r*r;
proc summary data=rDOpt;
   var psi2;
   weight w;
   output out=mPsi2 mean=mPsi2;
run;
title2 "Delta2";
proc means data=rDOpt max;
   var psi2;
run;
title2;
data _null_; set mPsi2;
   call symput('mPsi2',trim(left(mPsi2)));
data fig2004; set rDOpt(where=(w=0));
proc sort data=fig2004; by x;
run;

data dPoints; set dOpt;
   where (w > 0);
   dPsi2 = &mPsi2;
   rename x=dx;
data fig2004; merge fig2004 dPoints(keep=dx dPsi2);
run;

%let pagesize = %sysfunc(getoption(pagesize));
options ps=40;
proc plot data=fig2004;
   plot dPsi2*dx='*' psi2*x='.' / overlay;
run;
options ps=&pagesize;

title1;