/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 16.2                                              */
/*   TITLE: Table 16.2. Example 16.2: mixture experiment with linear  */
/*          constraints                                               */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Uses the ADX macros to create XVERT points.               */
/*          Also addresses SAS Tasks 16.4, 16.5, and 16.7.            */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Table 16.2. Example 16.2: mixture experiment with linear constraints";


/*
/  First approach: generate grid points and subset them.
/---------------------------------------------------------------------*/
data Grid; keep x1-x3;
   do ix1 = 0 to 100;
      do ix2 = 0 to 100;
         do ix3 = 0 to 100;
            if (ix1 + ix2 + ix3 = 100) then do;
               x1 = ix1 / 100;
               x2 = ix2 / 100;
               x3 = ix3 / 100;
               output;
               end;
            end; end; end;
data Grid; set Grid;
   where (  (0.2 <= x1 <= 0.7)
          & (0.1 <= x2 <= 0.6)
          & (0.1 <= x3 <= 0.6));
run;

/*
/  Find and print the optimum design.
/---------------------------------------------------------------------*/
proc optex data=Grid coding=orthcan;
   model x1|x2|x3@2 / noint;
   generate n=10 method=m_fedorov niter=1000 keep=10;
   output out=Design;
proc sort data=Design;
   by descending x1
      descending x2
      descending x3;
proc print data=Design;
run;



/*
/  Noting that all design points are vertices and edge cetroids, a
/  second approach uses XVERT to generate just these and to optimize
/  over them for different design sizes.  Begin by initializing the
/  ADX macros.
/---------------------------------------------------------------------*/
%adxgen;
%adxmix;

/*
/  Generate vertices and generalized edge centroids of the constrained
/  region, to be used as the candidate points.
/---------------------------------------------------------------------*/
%adxxvert(Candidates,x1 0.2-0.7 / x2 0.1-0.6 / x3 0.1-0.6);
data PointNumber; input PointNumber @@; datalines;
1 3 2 4 5 6 8 7 9 10 11 12 13
;
data Candidates; merge Candidates PointNumber;
run;

/*
/  A macro to find the optimum replications for the input candidates
/  for a given design size.  Save the replications and accumulate
/  them in a data set F for printing later.
/---------------------------------------------------------------------*/
%macro OptDesign(N);
proc optex data=Candidates coding=orthcan noprint;
   model x1|x2|x3@2 / noint;
   id PointNumber;
   generate n=&N method=m_fedorov niter=100 keep=10;
   output out=opDesign;
proc freq data=opDesign noprint;
   table PointNumber*x1*x2*x3 / out=f&N(rename=(Count=f&N) drop=PERCENT);
run;
%mend;
%OptDesign(10);
%OptDesign(20);
%OptDesign(30);
%OptDesign(35);
%OptDesign(100);
%OptDesign(10000);
data f; merge f10 f20 f30 f35 f100 f10000;
run;

/*
/  Compute the optimum weights over the candidates and add them to
/  the data set F as well.
/---------------------------------------------------------------------*/
proc iml;
   start DCritW(w) global(Z);
      Dw = diag((w##2)/sum(w##2));
      return(log(det(Z`*Dw*Z))/ncol(Z));
   finish;

   use Candidates; read all var ("x1":"x3") into xC;
   read all var ("PointNumber");
   Z = xC || xC[,1]#xC[,2] || xC[,1]#xC[,3] || xC[,2]#xC[,3];

   winit = j(nrow(Z),1) / nrow(Z);
   call nlpqn(rc,        /* Return code              */
              wOp,       /* Returned optimum weights */
              "DCritW",  /* Function to optimize     */
              winit,     /* Initial value of weights */
              1);        /* Specify a maximization   */
   wOp = (wOp##2)/sum(wOp##2);
   
   create w var {"PointNumber" "x1" "x2" "x3" "w"};
   data = PointNumber || xC || wOp`;
   append from data;

proc sort data=w; by PointNumber;
proc sort data=f; by PointNumber;
data f; merge f w; by PointNumber;
run;

proc print data=f noobs;
run;


/*
/  Erratum: The design given as optimum in Table 16.2 for N=35 is not
/  exactly so.
/---------------------------------------------------------------------*/
data ADT35; set Candidates;
   select (PointNumber);
      when ( 1) n = 4;
      when ( 2) n = 4;
      when ( 3) n = 3;
      when ( 4) n = 3;
      when ( 5) n = 3;
      when ( 6) n = 4;
      when ( 7) n = 0;
      when ( 8) n = 4;
      when ( 9) n = 4;
      when (10) n = 0;
      when (11) n = 4;
      when (12) n = 0;
      when (13) n = 2;
      end;
   do i = 1 to n; output; end;
run;

proc optex data=Candidates coding=orthcan;
   model x1|x2|x3@2 / noint;
   generate initdesign=ADT35 method=sequential;
   ods select Efficiencies;
   title2 "Evaluating the N=35 design in Table 16.2";
run;
   generate n=35 method=m_fedorov niter=100 keep=1;
   ods select Efficiencies;
   title2 "The true D-optimum N=35 design";
run;

title2;

title1;