/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 17.3                                              */
/*   TITLE: Fig. 17.7 & 17.8 - Sampling window for two consecutive    */
/*          first-order reactions                                     */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

%let True1 = 1.0;
%let True2 = 0.5;

title1 "Fig. 17.7 & 17.8 - Sampling window for two consecutive first-order reactions";

options nonotes;
proc iml;

   /*
   /  The Jacobian at given time points
   /------------------------------------------------------------------*/
   start f(Time);
      f1 = (  (&True2/((&True1-&True2)**2))
            + (&True1/ (&True1-&True2)    )*Time)#exp(-&True1*Time)
         - (  (&True2/((&True1-&True2)**2))     )#exp(-&True2*Time);
      f2 = (  (&True1/((&True1-&True2)**2))
            - (&True1/ (&True1-&True2)    )*Time)#exp(-&True2*Time)
         - (  (&True1/((&True1-&True2)**2))     )#exp(-&True1*Time);
      return(f1 || f2);
   finish;

   start DCrit(t);
      t1 = t[1]; t2 = t[2];
      F = f(t1) // f(t2);
      return((det((F`*F)/nrow(F)))**(1/ncol(F)));
   finish;

   /*
   /  Compute D-optimum 2-point exact design.
   /------------------------------------------------------------------*/
   tinit = {1 4};
   call nlpqn(rc,tOpt,"DCrit",tInit,1);
   tOpt = tOpt`;
   print tOpt (DCrit(tOpt));

   /*
   /  Standardized variance
   /------------------------------------------------------------------*/
   FOpt = f(tOpt);
   DOpt = DCrit(tOpt);
   start d(t) global(FOpt);
      f = f(t);
      return(f*inv(FOpt`*FOpt/nrow(FOpt))*f`);
   finish;

/*
/  d()/2 is the standardized variance relative to its maximum value,
/  which is p (=2 in this case).  diff() measures how far this is from
/  the target value of kappa.  
/---------------------------------------------------------------------*/
   start diff(t) global(kappa);
      return((d(t)/2 - kappa)**2);
   finish;

/*
/  The four points that define the sampling window are given by the
/  times above and below the optimum design times where d()/2 = kappa.
/  We use constrained nonlinear optimization to find these points for
/  given kappa=k and measure how far the minimum DCrit for these four
/  designs (and thus for the rectangle they define) is from the
/  target value of 0.9.
/---------------------------------------------------------------------*/
   start kCrit(k) global(tOpt,DOpt,kappa);
      kappa = k;
      call nlpqn(rc,t1L,"diff",tOpt[1] - 0.1,0,0       // tOpt[1]);
      call nlpqn(rc,t1U,"diff",tOpt[1] + 0.1,0,tOpt[1] // tOpt[:]);
      call nlpqn(rc,t2L,"diff",tOpt[2] - 0.1,0,tOpt[:] // tOpt[2]);
      call nlpqn(rc,t2U,"diff",tOpt[2] + 0.1,0,tOpt[2] // .      );
      return((min(   DCrit(t1L//t2L)/DOpt
                  || DCrit(t1U//t2L)/DOpt
                  || DCrit(t1L//t2U)/DOpt
                  || DCrit(t1U//t2U)/DOpt) - 0.9)**2);
   finish;
   call nlpqn(rc,k90,"kCrit",1);

/*
/  Resolve for given kappa=k90 to get the points, and print them and
/  their corresponding efficiencies.
/---------------------------------------------------------------------*/
   kappa = k90;
   call nlpqn(rc,t1L,"diff",tOpt[1] - 0.1,0,0       // tOpt[1]);
   call nlpqn(rc,t1U,"diff",tOpt[1] + 0.1,0,tOpt[1] // tOpt[:]);
   call nlpqn(rc,t2L,"diff",tOpt[2] - 0.1,0,tOpt[:] // tOpt[2]);
   call nlpqn(rc,t2U,"diff",tOpt[2] + 0.1,0,tOpt[2] // .      );
   tCorner = (t1L||t2L) // (t1U||t2L) // (t1L||t2U) // (t1U||t2U);
   DEff    = (   DCrit(tCorner[1,]`)
              // DCrit(tCorner[2,]`)
              // DCrit(tCorner[3,]`)
              // DCrit(tCorner[4,]`))/DOpt;
   print k90     [format=5.3] "  "
         tCorner [format=5.3] "  "
         DEff    [format=percent9.4];

title1;