/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 17.6                                              */
/*   TITLE: Example 17.6 and SAS Task 17.5 - Generalized consecutive  */
/*          reaction                                                  */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 17.6 and SAS Task 17.5 - A generalized consecutive reaction";

title2 "Table 17.4 - D-optimum designs for both rate and order";

data Region;
   do time = 0 to 20 by 0.01;
      output;
      end;
run;

%macro Table1704;

   data Table1704; if (0); run;
   %do Lambda2 = 1 %to 2;
   %do Lambda1 = 1 %to 2;

      proc model data=Region noprint;
         dependent C1 0.99999 C2 0.00001 C3 0 ;
         parms Theta1 1.0 Theta2 0.5 Lambda1 &Lambda1 Lambda2 &Lambda2;

/*
/ Equation 17.31 in the text: dCi/dt = Gi
/---------------------------------------------------------------------*/
         dert.C1 = -Theta1*(C1**Lambda1);
         dert.C2 =  Theta1*(C1**Lambda1) - Theta2*(C2**Lambda2);
         dert.C3 =  Theta2*(C2**Lambda2);

/*
/ dGi/dCj
/---------------------------------------------------------------------*/
         dg1C1 = -Theta1*Lambda1*(C1**(Lambda1-1));
         dg1C2 = 0;
         dg1C3 = 0;

         dg2C1 =  Theta1*Lambda1*(C1**(Lambda1-1));
         dg2C2 = -Theta2*Lambda2*(C2**(Lambda2-1));
         dg2C3 = 0;

         dg3C1 = 0;
         dg3C2 =  Theta2*Lambda2*(C2**(Lambda2-1));
         dg3C3 = 0;

/*
/ dGi/dPj
/---------------------------------------------------------------------*/
         dg1p1 = -C1**(Lambda1);
         dg1p2 = 0;
         dg1p3 = -Theta1*(C1**Lambda1)*log(C1);
         dg1p4 = 0;

         dg2p1 =  C1**Lambda1;
         dg2p2 = -C2**Lambda2;
         dg2p3 =  Theta1*(C1**Lambda1)*log(C1);
         dg2p4 = -Theta2*(C2**Lambda2)*log(C2);

         dg3p1 = 0;
         dg3p2 =  C2**Lambda2;
         dg3p3 = 0;
         dg3p4 =  Theta2*(C2**Lambda2)*log(C2);
      
/*
/ Equation 17.36 in the text.
/---------------------------------------------------------------------*/
         dert.f11 = dg1C1*f11 + dg1C2*f21 + dg1C3*f31 + dg1p1;
         dert.f21 = dg2C1*f11 + dg2C2*f21 + dg2C3*f31 + dg2p1;
         dert.f31 = dg3C1*f11 + dg3C2*f21 + dg3C3*f31 + dg3p1;

         dert.f12 = dg1C1*f12 + dg1C2*f22 + dg1C3*f32 + dg1p2;
         dert.f22 = dg2C1*f12 + dg2C2*f22 + dg2C3*f32 + dg2p2;
         dert.f32 = dg3C1*f12 + dg3C2*f22 + dg3C3*f32 + dg3p2;

         dert.f13 = dg1C1*f13 + dg1C2*f23 + dg1C3*f33 + dg1p3;
         dert.f23 = dg2C1*f13 + dg2C2*f23 + dg2C3*f33 + dg2p3;
         dert.f33 = dg3C1*f13 + dg3C2*f23 + dg3C3*f33 + dg3p3;

         dert.f14 = dg1C1*f14 + dg1C2*f24 + dg1C3*f34 + dg1p4;
         dert.f24 = dg2C1*f14 + dg2C2*f24 + dg2C3*f34 + dg2p4;
         dert.f34 = dg3C1*f14 + dg3C2*f24 + dg3C3*f34 + dg3p4;

         solve C1 C2 C3 f11 f21 f31
                        f12 f22 f32
                        f13 f23 f33
                        f14 f24 f34 / out=Candidates time=time;
      run;

      proc optex data=Candidates coding=orthcan noprint;
         model f21 f22 f23 f24 / noint;
         generate n=4 method=m_fedorov niter=1000 keep=10;
         id Time;
         output out=Design;
      proc sort data=Design; by time;
      proc transpose data=Design out=tDesign(rename=(COL1=t1 COL2=t2 COL3=t3 COL4=t4));
         var time;
      data tDesign; set tDesign;
         Lambda1 = &Lambda1;
         Lambda2 = &Lambda2;
      data Table1704; set Table1704 tDesign;
      run;
      %end; %end;

   proc print data=Table1704 noobs;
      var Lambda1 Lambda2 t:;
   run;
%mend;
%Table1704;

title2;

title1;