/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 17.7                                              */
/*   TITLE: Example 17.7 - Multivariate nonlinear design              */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Uses NLIN to create the derivatives as callable code, in  */
/*          order to find optimum design support and weights using    */
/*          only the input models.                                    */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 17.7 - Multivariate nonlinear design";

/*
/  Models for four chemical concentrations; see equations (17.41) in
/  the text.
/---------------------------------------------------------------------*/
%let ModelA = exp(-t1*t);
%let ModelB = 
        ((     t3     /( t2+t3            )) * (1          - exp(-(t2+t3)*t)))
      - ((((t1-t3)*t4)/((t1-t2-t3)*(t4-t1))) * (exp(-t1*t) - exp(-(t2+t3)*t)))
      - ((((t4-t3)*t1)/((t2+t3-t4)*(t4-t1))) * (exp(-t4*t) - exp(-(t2+t3)*t)));
%let ModelD = (t4/(t4-t1)) * (exp(-t1*t) - exp(-t4*t));
%let ModelC = 1 - (&ModelA) - (&ModelB) - (&ModelD);


/*
/  The following code creates macro variables with SAS statements that
/  compute the derivatives, as variables named D<<Response>>_t<<i>>,
/  for given <<Response>> and for <<i>> = 1,2,3,4.  First the LIST
/  option in PROC NLIN creates a data set Deriv with the statements
/  for each response's derivatives
/---------------------------------------------------------------------*/
%macro MakeDeriv(Response,Model);
   %global D1 D2 D3 D4;

   data _temp; &Response = 1; t = 1; run;
   ods listing close;
   proc nlin data=_temp list;
      parameters t1=1 t2=1 t3=1 t4=1;
      model &Response = &Model;
      ods output ProgList=Deriv;
   run;
   ods listing;
   data Deriv; set Deriv;
      rx = rxparse("'MODEL.' to ''");
      call rxchange(rx,999,Source);
      rx = rxparse("'/@' to '_D'");
      call rxchange(rx,999,Source);
      rx = rxparse("'@' to 'D'");
      call rxchange(rx,999,Source);
      drop rx;
   run;
%mend;

data Code; if (0); run;
%MakeDeriv(A,&ModelA); data Code; set Code Deriv; run;
%MakeDeriv(B,&ModelB); data Code; set Code Deriv; run;
%MakeDeriv(C,&ModelC); data Code; set Code Deriv; run;
%MakeDeriv(D,&ModelD); data Code; set Code Deriv; run;

%let nCode = 0;
data _null_; set Code;
   call symput('nCode',trim(left(symget('nCode')+1)));
   call symput('Code'||trim(left(symget('nCode'))),Source);
run;



/*
/  Now, create a macro program that uses these code macros to define
/  an IML function that computes the sensitivities for a given single
/  time value.
/---------------------------------------------------------------------*/
%macro DefineMake1Z;
start Make1Z(t) global(theta);
   t1 = theta[1];
   t2 = theta[2];
   t3 = theta[3];
   t4 = theta[4];
   DA_Dt1 = 0; DA_Dt2 = 0; DA_Dt3 = 0; DA_Dt4 = 0;
   DB_Dt1 = 0; DB_Dt2 = 0; DB_Dt3 = 0; DB_Dt4 = 0;
   DC_Dt1 = 0; DC_Dt2 = 0; DC_Dt3 = 0; DC_Dt4 = 0;
   DD_Dt1 = 0; DD_Dt2 = 0; DD_Dt3 = 0; DD_Dt4 = 0;
   %do iCode = 1 %to &nCode;
      &&Code&iCode
      %end;
   XA = DA_Dt1 || DA_Dt2 || DA_Dt3 || DA_Dt4;
   XB = DB_Dt1 || DB_Dt2 || DB_Dt3 || DB_Dt4;
   XC = DC_Dt1 || DC_Dt2 || DC_Dt3 || DC_Dt4;
   XD = DD_Dt1 || DD_Dt2 || DD_Dt3 || DD_Dt4;
   return(XA || XB || XC || XD);
finish;
%mend;




/*
/  The schema for finding optimum weights and support points is much
/  the same as we have used in many previous programs: find optimal
/  weights over a grid, then jointly optimize weights and support.
/  The new wrinkle is that the DCrit() function depends on which
/  reponses are assumed to be observed.
/---------------------------------------------------------------------*/
%let Truet1 = 0.8;
%let Truet2 = 0.2;
%let Truet3 = 0.15;
%let Truet4 = 1.4;
proc iml;

   %DefineMake1Z;
   start MakeZ(x) global(theta);
      t = shape(x,nrow(x)*ncol(x),1);
      free Z;
      do i = 1 to nrow(t);
         Z = Z // Make1Z(t[i]);
         end;
      return(Z);
   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,nrow(wx)*ncol(wx)/2)[1,];
      w = (w##2)/sum(w##2);
      x = shape(wx,2,nrow(wx)*ncol(wx)/2)[2,];
      return(DCrit(w`,x`));
   finish;

   start DCrit(_w,x) global(which);
      XAll = MakeZ(x);
      X1 = XAll[,0*4 + (1:4)];
      X2 = XAll[,1*4 + (1:4)];
      X3 = XAll[,2*4 + (1:4)];
      X4 = XAll[,3*4 + (1:4)];
      W = diag(_w);
      M = 1e-12*i(4);
      if (index(which,"A")) then M = M + X1`*W*X1;
      if (index(which,"B")) then M = M + X2`*W*X2;
      if (index(which,"C")) then M = M + X3`*W*X3;
      if (index(which,"D")) then M = M + X4`*W*X4;
      dd = llog(det(M));
      return(dd);
   finish;

   start dd(wD,xD,x) global(which);
      XAll = MakeZ(xD);
      X1 = XAll[,0*4 + (1:4)];
      X2 = XAll[,1*4 + (1:4)];
      X3 = XAll[,2*4 + (1:4)];
      X4 = XAll[,3*4 + (1:4)];
      W = diag(wD);
      M = 0*i(4);
      if (index(which,"A")) then M = M + X1`*W*X1;
      if (index(which,"B")) then M = M + X2`*W*X2;
      if (index(which,"C")) then M = M + X3`*W*X3;
      if (index(which,"D")) then M = M + X4`*W*X4;

      F = MakeZ(x);
      f1 = F[,0*4 + (1:4)];
      f2 = F[,1*4 + (1:4)];
      f3 = F[,2*4 + (1:4)];
      f4 = F[,3*4 + (1:4)];
      dd = 0;
      if (index(which,"A")) then dd = dd + f1*ginv(M)*f1`;
      if (index(which,"B")) then dd = dd + f2*ginv(M)*f2`;
      if (index(which,"C")) then dd = dd + f3*ginv(M)*f3`;
      if (index(which,"D")) then dd = dd + f4*ginv(M)*f4`;
      return(dd);
   finish;

   theta = {&Truet1 &Truet2 &Truet3 &Truet4};

%macro OptDesign(which);
   which = &which;

   nGrid = 21;
   x    = 0 + (20 - 0)*((1:nGrid)`-1)/(nGrid-1);
   winit = j(nrow(x),1) / nrow(x);

   call nlpqn(rc,wGrid,"DCritW",winit,1);
   wGrid = ((wGrid##2)/sum(wGrid##2))`;
   xGrid = x  [loc(wGrid > 1e-3)]`;
   wGrid = wGrid[loc(wGrid > 1e-3)]`; wGrid = wGrid /sum(wGrid);

   xOp = xGrid;
   wOp = wGrid;

   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(wOp)) || shape( 0,1,ncol(xOp)))
            // (shape(.,1,ncol(wOp)) || shape(20,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-3,1e-4);

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


   /*
   /  Check the design for optimality by seeing whether the variance
   /  is maximized at design points, and is equal to p.
   /------------------------------------------------------------------*/
   free tdd;
   do i = 1 to ncol(xOp);
      tdd = tdd // (xOp[i] || dd(wOp`,xOp`,xOp[i]));
      end;
   maxt  = .;
   maxdd = -1e24;
   do t = 0 to 20 by 0.1;
      if (min(abs(xOp - t)) > 1e-2) then do;
         dd = dd(wOp`,xOp`,t);
         if (dd > maxdd) then do;
            maxt  = t;
            maxdd = dd;
            end;
         end;
      end;

   /*
   /  Discover the number of parameters p by looking at the non-zero
   /  sensitivities at an arbitrary time (t=2, in this case).
   /------------------------------------------------------------------*/
   X = abs(MakeZ(1.2345));
   free pnz;
   if (index(which,"A")) then pnz = pnz // X[<>,0*4 + (1:4)];
   if (index(which,"B")) then pnz = pnz // X[<>,1*4 + (1:4)];
   if (index(which,"C")) then pnz = pnz // X[<>,2*4 + (1:4)];
   if (index(which,"D")) then pnz = pnz // X[<>,3*4 + (1:4)];
   p = sum(pnz[<>,] > 0);
   
   Check = ((max(abs(p - tdd[,2])) < 1e-4) & (maxdd <= min(tdd[,2])));

   /*
   /  Also find the optimum four-point design, to compare to the
   /  D-optimum design.
   /------------------------------------------------------------------*/
   nD = 4;
   f = round(4*wOp);
   if (sum(f) = nD) then do;
      free xinit;
      do i = 1 to ncol(wOp);
         do j = 1 to f[i];
            xinit = xinit || xOp[i];
            end;
         end;
      end;
   else do;
      xinit = 0.1 + (19.9 - 0.1)*((1:nD)`-1)/(nD-1);
      end;
   w     = j(nD,1)/nD;
   con = shape( 0,1,nD) // shape(20,1,nD);
   call nlpqn(rc,xEq,"DCritX",xinit,1,con);
   xEq = xEq`;
   if (inputn(symget('sysver'),'3.1') >= 9) then do;
      start sort(x,idx);
         sortx = j(nrow(x), ncol(x), .);
         do i = nrow(idx)*ncol(idx) to 1 by -1;
            sortx[rank(x[,idx[i]]),] = x[];
            end;
         x = sortx;
      finish;
      end;
   call sort(xEq,1);
   bestD   = DCritX(xEq);
   bestxEq = xEq;
   do i = 1 to 100;
      xinit = 0 + (20 - 0)*ranuni(1:nD);
      call nlpqn(rc,xEq,"DCritX",xinit,1,con);
      D = DCritX(xEq);
      xEq = xEq`;
      call sort(xEq,1);
      if (D > bestD) then do;
         bestD = D;
         bestxEq = xEq;
         end;
      end;

   Weight  = wOp`;
   Support = xOp`;
   EqEff = 100*exp(DCrit(w   ,BestxEq)/4)
              /exp(DCrit(wOp`,xOp`   )/4);
   print "Optimum design for (" &which ")",,
      Weight  [format=8.4]
      Support [format=8.4]
      BestxEq [format=8.4]
      EqEff   [format=8.4];
   if (^Check) then do;
      ddDesign = tdd[,2];
      print "Problem: " ddDesign Maxdd;
      end;
   else do;
      end;
%mend;

options nonotes;
%OptDesign("A   ");
%OptDesign(" B  ");
%OptDesign("  C ");
%OptDesign("   D");
%OptDesign("AB  ");
%OptDesign("A C ");
%OptDesign("A  D"); /* Something wrong with this one */
%OptDesign(" BC ");
%OptDesign(" B D");
%OptDesign("  CD");
%OptDesign("ABC ");
%OptDesign("AB D");
%OptDesign("A CD");
%OptDesign(" BCD");
%OptDesign("ABCD");

title1;