/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 24.2                                              */
/*   TITLE: Example 24.2 - Compartmental model                        */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Not quite the designs given in the text, but equivalent   */
/*          with respect to the optimality criterion, as demonstrated */
/*          by the printed efficiency comparisons.                    */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 24.2 - Compartmental model";
title2 "Table 24.1 - Optimum multiple series designs for theophylline sampling";

%let Truet1  =  4.2899;
%let Truet2  =  0.05885;
%let Truet3  = 21.7965;
%let TrueSig = 0.9988;
options nonotes;
proc iml;

   start y(Time) global(t1,t2,t3);
      B = t3*(exp(-t1*Time) - exp(-t2*Time));
      return(B);
   finish;

   start MakeZ(Time) global(t1,t2,t3);
      Z =     t3*                      Time#exp(-t1*Time)
          || -t3* Time#exp(-t2*Time)
          ||     (     exp(-t2*Time) -      exp(-t1*Time));
      return(Z);
   finish;
   
   start MakeW(Time) global(rho);
      n = nrow(Time)*ncol(Time);
      if (rho = .) then do;
         return(i(n));
         end;
      else do;
         W = shape(.,n,n);
         do i = 1 to n;
            do j = 1 to n;
              W[i,j] = exp(-rho*abs(Time[i]-Time[j]));
               end;
            end;
         end;
      return(W);
   finish;

   start llog(x); if (x > 0) then return(log(x)); else return(-9999); finish;

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

   start M1(_x);
      n = nrow(_x)*ncol(_x);
      x = shape(_x,n,1);
      F = MakeZ(x);
      W = MakeW(x);
      return(F`*inv(W)*F/nrow(F));
   finish;

   start DCrit(w) global(Des);
      M = w[1]*M1(Des[1,]);
      do i = 2 to nrow(w)*ncol(w);
         M = M + w[i]*M1(Des[i,]);
         end;
      return(llog(det(M))/ncol(M));
   finish;

   t1  = &Truet1;
   t2  = &Truet2;
   t3  = &Truet3;

   start dDCrit(w,x);
      F = MakeZ(x);
      return(det(F`*diag(w)*F)**(1/ncol(F)));
   finish;

   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;

   %macro MakeDesign(rho,OuterInc,nInnerInc,InnerInc,nIter=1,print=1);
      rho = ρ

      free Can;
      do t = 0 to 20-(&nInnerInc-1)*&InnerInc by &OuterInc;
         free Cani;
         do j = 1 to &nInnerInc;
            Cani = Cani || t+(j-1)*&InnerInc;
            end;
       if ((Cani[1] < 3) | (Cani[&nInnerInc] > 16)) then  
            Can = Can // Cani;
         end;

      free All;
      BestD = -9999;
      do iInit = 0 to &nIter;
         Des = Can;
         if (iInit = 0) then wOp = j(1,nrow(Des))/nRow(Des);
         else                wOp = ranuni(1:nrow(Des))`;
         Changed = 1;
         iter = 0;
         do while (Changed & (iter < 100));
            wLast = wOp;
            winit = sqrt(wOp);
            call nlpqn(rc,wOp,"DCritW",winit,1);
            wOp = ((wOp##2)/sum(wOp##2))`;
            Des = Des[loc(wOp > 1e-3),];
            wOp = wOp[loc(wOp > 1e-3) ]; wOp = wOp /sum(wOp);
            Changed = 0;
            if      (nrow(wOp) ^= nrow(wLast))     then Changed = 1;
            else if (max(abs(wOp - wLast)) > 1e-8) then Changed = 1;
            iter = iter + 1;
            end;
         d = DCrit(wOp);
         if (d > BestD) then do;
            BestD   = d;
            BestW   = wOp;
            BestDes = Des;
            end;
         All = All // (iInit || d);
         end;
      wOp = BestW;
      xOp = BestDes;
      All[,2] = -All[,2];
      call sort(All,2);
      %if (&print >= 2) %then %do;
         print All;
         %end;
      All[,2] = -All[,2];

      %if (&print >= 1) %then %do;
         xCmp = xCmp1` || xCmp1`+0.25;
         wCmp = wCmp` / sum(wCmp);
*         print wOp xOp "    " wCmp xCmp;
         Des = xOp;  dOp  = DCrit(wOp );
         Des = xCmp; dCmp = DCrit(wCmp);
         tOp = substr(putn(60*xOp,'time.'),4);
         print ("Tau = "+trim(left(char(rho,12)))),,
               (tOp`)              ({.})                    [format=8.4],,
               (wOp`) [format=5.3] (100*exp(dCmp)/exp(dOp)) [format=8.4];
         %end;
   %mend;

%let Time0 = %sysfunc(datetime());

xCmp1 = {0     0.25  17.75 };
wCmp  = {0.239 0.428  0.333};
%MakeDesign( 0.5,0.25,2,0.25,nIter=0,print=1);

xCmp1 = {0     1.25  18.25 };
wCmp  = {0.333 0.333  0.333};
%MakeDesign( 1.0,0.25,2,0.25,nIter=0,print=1);

xCmp1 = {0.25  1.25  18.25 };
wCmp  = {0.400 0.267  0.333};
%MakeDesign( 4.0,0.25,2,0.25,nIter=0,print=1);

xCmp1 = {0.25  1.25   1.5   18.25};
wCmp  = {0.373 0.106  0.188 0.333};
%MakeDesign( 8.0,0.25,2,0.25,nIter=0,print=1);

xCmp1 = {0.25  1.5   18.5  };
wCmp  = {0.364 0.303  0.333};
%MakeDesign(16.0,0.25,2,0.25,nIter=0,print=1);

%let Time1 = %sysfunc(datetime());
%put Time = %sysevalf(&Time1 - &Time0);

title2;
title1;