/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 16.1                                              */
/*   TITLE: Figures 16.1 & 16.2 - Optimal mixture designs             */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Uses the ADX macros to create simplex lattice designs.    */
/*          Also addresses SAS Task 16.3.                             */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Figures 16.1 & 16.2 - Optimal mixture designs";

/*
/  Initiate ADX macros.
/---------------------------------------------------------------------*/
%adxgen;
%adxmix;

/*
/  The default order for a simplex lattice design is two.
/---------------------------------------------------------------------*/
%adxsld(SLD2,x1 x2 x3);

/*
/  The third parameter of ADXSLD changes the order.
/---------------------------------------------------------------------*/
%adxsld(SLD3,x1 x2 x3,3);

/*
/  Use OPTEX with a standard cubic model in two of the mixture factors
/  to approximate the exact optimum design on a large simplex lattice
/  grid.
/---------------------------------------------------------------------*/
%adxsld(Grid,x1 x2 x3,100);
proc optex data=Grid coding=orthcan;
   title2 "Figure 16.2a - Efficiencies for third-order simplex lattice design";
   model x1 x2 x1*x1 x1*x2 x2*x2 x1*x1*x1 x1*x1*x2 x1*x2*x2 x2*x2*x2;
   generate initdesign=SLD3 method=sequential;
run;
   title2 "Figure 16.2b1 - Optimal cubic design over order-100 simplex lattice grid";
   generate n=10 method=m_fedorov niter=100 keep=1;
   id x3;
   output out=opDesign;
proc sort data=opDesign;
   by x1 x2 x3;
run;
title2;



/*
/  Use the nonlinear optimization facilities in IML to optimally
/  improve the SLD3.
/---------------------------------------------------------------------*/
proc iml;

   start MakeZ(xx);
      x = xx || 1 - xx[,+];
      Z = x;
      do i = 1 to ncol(x)-1; do j = i+1 to ncol(x);
         Z = Z || x[,i]#x[,j];
         end; end;
      do i = 1 to ncol(x)-1; do j = i+1 to ncol(x);
         Z = Z || x[,i]#x[,j]#(x[,i]-x[,j]);
         end; end;
      do i = 1 to ncol(x)-2; do j = i+1 to ncol(x)-1; do k = j+1 to ncol(x);
         Z = Z || x[,i]#x[,j]#x[,k];
         end; end; end;
      return(Z);
   finish;

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

   start DCritX(x);
      xx = shape(x,nrow(x)*ncol(x)/2,2);
      Z = MakeZ(xx);
      return(llog(det(Z`*Z))/ncol(Z));
   finish;

   /*
   /  The hard part about this approach is setting up the right
   /  constraints.  We optimize on the first two mixture components
   /  only, to obviate the equality constraint.  So we need to set up
   /  lower and upper constraints not only in x1 and x2, but also on
   /  their sum.
   /------------------------------------------------------------------*/
   use SLD3;
   read all var ("x1":"x3") into xSLD3;
   Z = MakeZ(xSLD3[,1:2]);
   N = nrow(xSLD3);
   xinit = shape(xSLD3[,1:2],2*N,1);
   con =    (shape( 0,1,2*N) ||       { . .}     )    /* xij     >= 0 */
         // (shape( 1,1,2*N) ||       { . .}     )    /* xij     <= 1 */
         // ((i(N)@j(1,2))   || shape({ 1 0},N,2))    /* xi1+xi2 >= 0 */
         // ((i(N)@j(1,2))   || shape({-1 1},N,2));   /* xi1+xi2 <= 1 */
   call nlpqn(rc,xx,"DCritX",xinit,1,con);
   xim = shape(xx,nrow(xx)*ncol(xx)/2,2);
   xim = xim || 1 - xim[,+];
   create Improved var {"x1" "x2" "x3"}; append from xim;

   /*
   /  Although we optimized the support assuming equal weights, this
   /  is in fact the optimum design.  We demonstrate this by
   /  evaluating the prediction variance over the large grid and
   /  showing that the maximum is 10, the number of parameters.
   /------------------------------------------------------------------*/
   use opDesign;
   read all var ("x1":"x3") into xop;

   Z = MakeZ(xSLD3[,1:2]); VSC = inv(Z`*Z/nrow(Z));
   Z = MakeZ(xim  [,1:2]); Vim = inv(Z`*Z/nrow(Z));
   Z = MakeZ(xop  [,1:2]); Vop = inv(Z`*Z/nrow(Z));
   use Grid;
   read all var ("x1":"x3") into xCan;
   Z = MakeZ(xCan[,1:2]);
   free var;
   do i = 1 to nrow(Z);
      var = var // (Z[i,]*VSC*Z[i,]` || Z[i,]*Vop*Z[i,]` || Z[i,]*Vim*Z[i,]`);
      end;
   GEff = ncol(Z)/var[<>,];

   /*
   /  Compute the D-efficiencies and OPTEX's form of the G-efficiency,
   /  too: recall that OPTEX defines the G-efficiency as the square
   /  root of the usual form.  Compare the OPTEX G-efficiencies to the
   /  ones computed by OPTEX.
   /------------------------------------------------------------------*/
   DEff =    exp(DCritX(xSLD3[,1:2]))
          || exp(DCritX(xOp  [,1:2]))
          || exp(DCritX(xIm  [,1:2]));
   DEff = DEff / max(DEff);

   print "Design Efficiencies",,
         (DEff // GEff // sqrt(GEff)) [format=percent10.4
                         rowname={"D-Efficiency" "G-Efficiency" "OPTEX G-Eff"}
                         colname={"SLD3" "OPTEX" "D-Optimum"}];



title2 "Figure 16.1 - Second-order simplex lattice design";
proc print data=SLD2;
run;
title2;

title2 "Figure 16.2a - Third-order simplex lattice design";
proc print data=SLD3;
run;
title2;

title2 "Figure 16.2b1 - Optimal cubic design over order-100 simplex lattice grid";
proc print data=opDesign;
run;
title2;

title2 "Figure 16.2b2 - Optimum cubic design";
proc print data=Improved;
run;
title2;

title1;