/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 18.5                                              */
/*   TITLE: Tables 18.12 - Local and Bayesian D-optimum designs for   */
/*          reversible reaction                                       */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: The flatness of this optimization problem makes it very   */
/*          sensitive numerically.  Different machines will get       */
/*          different answers, and even on the same machine, slight   */
/*          changes in                                                */
/*             * how the computations are performed,                  */
/*             * how the optimization is initialized, and             */
/*             * what values are used for consolidation               */
/*          can change the results.  It took ~2.5 minutes to run on   */
/*          server-class PC in May, 2007.                             */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 18.4 - Reversible reaction";

proc iml;

/*
/  For this problem, instead of simply using the direct code for the
/  derivatives, we do our best to vectorize and optimize the
/  calculation.
/---------------------------------------------------------------------*/
   start MakeZ(x,theta);
      t = shape(x,nrow(x)*ncol(x),1);
      t1 = theta[1];
      t2 = theta[2];
      t3 = theta[3];
      et1t       = exp(-t1*t);
      et2t3      = exp(-(t2+t3)*t);
      tet1t      = t#et1t;
      tet2t3     = t#et2t3;
      t1t2t3     = t1-t2-t3;
      t1t3       = t1-t3;
      sqt1t2t3   = t1t2t3**2;
      sqt2t3     = (t2+t3)**2;
      et1tmet2t3 = et1t-et2t3;
      onemet2t3  = 1.0-et2t3;

      temp1 = t1t3/sqt1t2t3*et1tmet2t3;
      temp2 =   t3/(t2+t3) *tet2t3;
      temp3 =   t3/sqt2t3  *onemet2t3;
      temp4 =    1/t1t2t3  *et1tmet2t3;
      temp5 = t1t3/t1t2t3  *tet2t3;

      X2 =    -temp4 + temp1 + t1t3/t1t2t3*tet1t
           || -temp3 + temp2 - temp1 - temp5
           || 1/(t2+t3)*onemet2t3 - temp3 + temp2 + temp4 - temp1 - temp5;

   return(X2);
   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);
      if (any(w < 1e-4)) then w[loc(w < 1e-4)] = 0;
      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(STheta);
      dd = 0;
      do iTheta = 1 to nrow(STheta);
         X2 = MakeZ(x,STheta[iTheta,]);
         dd = dd + llog(det(1e-8*i(3) + X2`*diag(_w)*X2));
         end;
      return(dd/nrow(sTheta));
   finish;

/*
/  Generate a random sample of parameter values of size nTheta from
/  a multivariate lognormal distribution with the given variance tau.
/---------------------------------------------------------------------*/
%macro MakeSample(tau,nTheta);
   %let Truet1 = 0.7;
   %let Truet2 = 0.2;
   %let Truet3 = 0.15;
   tTheta = {&Truet1 &Truet2 &Truet3};
   tau = &tau;
   nTheta = &nTheta;
   free Stheta;
   do while (nrow(sTheta) < nTheta);
      theta = exp(log(tTheta) + tau*rannor(1:3));
      if (1.1*(theta[2] + theta[3]) <= theta[1]) then
         Stheta = Stheta // theta;
      end;
%mend;

/*
/  Find the D-optimum design for the current sample of parameter
/  values.
/---------------------------------------------------------------------*/
%macro Optimize;
   nint = 21;
   logx  = log(0.01) + (log(20) - log(0.01))*((1:nint)`-1)/(nint-1);
   x     = exp(logx);
   winit = j(nrow(x),1) / nrow(x);
   call nlpqn(rc,wOp,"DCritW",winit,1);
   wOp = ((wOp##2)/sum(wOp##2));
   xOp = x;
   call Consolidate(xOp,wOp,0,1e-3);
   lastDCrit = DCrit(wOp`,xOp`);

   i = 0;
   diff = 1;
   do while ((i < 100) & (diff > 1e-8));
      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-1,5e-3);
      thisDCrit = DCrit(wOp`,xOp`);

      if (ncol(xOp) ^= nrow(last)) then
         diff = 1;
      else
         diff = abs(lastDCrit - thisDCrit);
      lastDCrit = thisDCrit;
      end;
%mend;

options nonotes;

%MakeSample(0.5*log(4), 50); %Optimize; wOp1 = wOp`; xOp1 = xOp`;
%MakeSample(0.5*log(4), 50); %Optimize; wOp2 = wOp`; xOp2 = xOp`;
%MakeSample(0.5*log(4), 50); %Optimize; wOp3 = wOp`; xOp3 = xOp`;
%MakeSample(0.5*log(4),500); %Optimize; wOp0 = wOp`; xOp0 = xOp`;
%MakeSample(0.5*log(1),  1); %Optimize; wOp4 = wOp`; xOp4 = xOp`;
%MakeSample(0.5*log(2), 50); %Optimize; wOp5 = wOp`; xOp5 = xOp`;

wx = shape(.,6*2,nrow(xOp0));
wx[ 1: 2,1:(nrow(xOp4)-1)] = (xOp4` // wOp4`)[,1:(nrow(xOp4)-1)];
wx[ 1: 2,   nrow(xOp0)   ] = (xOp4` // wOp4`)[,   nrow(xOp4)   ];
wx[ 3: 4,1:(nrow(xOp5)-1)] = (xOp5` // wOp5`)[,1:(nrow(xOp5)-1)];
wx[ 3: 4,   nrow(xOp0)   ] = (xOp5` // wOp5`)[,   nrow(xOp5)   ];
wx[ 5: 6,1:(nrow(xOp1)-1)] = (xOp1` // wOp1`)[,1:(nrow(xOp1)-1)];
wx[ 5: 6,   nrow(xOp0)   ] = (xOp1` // wOp1`)[,   nrow(xOp1)   ];
wx[ 7: 8,1:(nrow(xOp2)-1)] = (xOp2` // wOp2`)[,1:(nrow(xOp2)-1)];
wx[ 7: 8,   nrow(xOp0)   ] = (xOp2` // wOp2`)[,   nrow(xOp2)   ];
wx[ 9:10,1:(nrow(xOp3)-1)] = (xOp3` // wOp3`)[,1:(nrow(xOp3)-1)];
wx[ 9:10,   nrow(xOp0)   ] = (xOp3` // wOp3`)[,   nrow(xOp3)   ];
wx[11:12,1:(nrow(xOp0)-1)] = (xOp0` // wOp0`)[,1:(nrow(xOp0)-1)];
wx[11:12,   nrow(xOp0)   ] = (xOp0` // wOp0`)[,   nrow(xOp0)   ];

tau = {"0.5*log(1)  "
       "0.5*log(2)  "
       "0.5*log(4) 1"
       "0.5*log(4) 2"
       "0.5*log(4) 3"
       "0.5*log(4)  "};
nTheta = {1 50 50 50 50 500};

   print "Table 18.12: Locally and Bayesian D-Optimal Designs",,
      "   tau  " "    n " "---------------------- t // w ----------------------",,
      (tau[1]) (nTheta[1]) [format=best3.] (wx[ 1: 2,]) [format=8.4],,
      (tau[2]) (nTheta[2]) [format=best3.] (wx[ 3: 4,]) [format=8.4],,
      (tau[3]) (nTheta[3]) [format=best3.] (wx[ 5: 6,]) [format=8.4],,
      (tau[4]) (nTheta[4]) [format=best3.] (wx[ 7: 8,]) [format=8.4],,
      (tau[5]) (nTheta[5]) [format=best3.] (wx[ 9:10,]) [format=8.4],,
      (tau[6]) (nTheta[6]) [format=best3.] (wx[11:12,]) [format=8.4];

title1;