/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 20.2                                              */
/*   TITLE: Example 2.2 & SAS Task 20.2 - Second-order response       */
/*          surface: constrained design region                        */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: This program generates optimum continuous designs by      */
/*          optimizing only weights for given support, comparing      */
/*          results to optimum replications for large exact designs   */
/*          over the same support.                                    */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 20.2 - Second-order response surface: constrained design region";

data Grid;
   do ix1 = 1 to 21;
      do ix2 = 1 to 21;
         x1 = -1 + 2*(ix1-1)/20;
         x2 = -1 + 2*(ix2-1)/20;
         if (  (2*x1 + x2 <=  1)
             & (  x1 + x2 >= -1)
             & (x2 - x1 <=  1.5)) then output;
         end;
      end;
run;

%macro N0(alpha,N);                       /* Inverting equation 20.5: */
   %if (%sysevalf(&alpha > 0)) %then %do; /* prior N0 corresponding   */
      %sysevalf(((1-&alpha)/&alpha)*&N)   /* to given alpha           */
      %end;
   %else %do;
      %sysevalf(1000*&N) 
      %end;
%mend;

/*
/  Use OPTEX to generate large exact designs for given alpha.
/---------------------------------------------------------------------*/
%macro Exact(N,alpha,ds,wgt);
   proc optex data=Grid coding=none noprint;
      model x1 x2, x1*x1 x1*x2 x2*x2 / prior = 0,%N0((&alpha),&N);
      generate n=&N niter=10 keep=10 method=m_fedorov;
      output out=Design;
   proc freq data=Design noprint;
      table x2*x1 / list nocum out=&ds(rename=(Percent=&wgt));
   run;
%mend;
%Exact(1000,0.318,f318,f318);
%Exact(1000,0.531,f531,f531);
%Exact(1000,0.647,f647,f647);
%Exact(1000,0.949,f949,f949);



/*
/  Optimal weights for given alpha.
/---------------------------------------------------------------------*/
proc iml;
   start MakeZ(_x);
      x = shape(_x,ncol(_x)*nrow(_x)/2,2);
      Z = j(nrow(x),1) || x[,1] || x[,2]
                       || x[,1]#x[,1]
                       || x[,1]#x[,2]
                       || x[,2]#x[,2];
      R12 = I(nrow(Z)) - Z[,1:3]*ginv(Z[,1:3]`*Z[,1:3])*Z[,1:3]`;
      Z = Z[,1:3] || R12*Z[,4:6];
      return(Z);
   finish;

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

   start DCritW(w) global(Z,alpha);
      d = det(   (1-alpha)   *diag((1:ncol(Z))>=4)
               +    alpha *Z`*diag((w##2)/sum(w##2))*Z);
      return(llog(d));
   finish;

/*
/  Quick and dirty, knowing that the support of the optimum design is
/  contained in the candidate set.
/---------------------------------------------------------------------*/
   start AlphaOpt;
      x     = xCan; Z = MakeZ(x);            /* Optimize weights for  */
      winit = sqrt(j(nrow(x),1) / nrow(x));  /* all support points.   */
      call nlpqn(rc,w,"DCritW",winit,1);
      w = (w##2)`/sum(w##2);
      xOp = x[loc(w > 1e-5),];               /* Trim support points   */
      wOp = w[loc(w > 1e-5) ];               /* with negligible wgts. */
      wOp = wOp /sum(wOp);

      x = xOp; Z = MakeZ(x);                 /* And optimize again.   */
      winit = wOp;
      call nlpqn(rc,w,"DCritW",winit,1);
      w = ((w##2)/sum(w##2))`;
      xOp = x[loc(w > 1e-5),];
      wOp = w[loc(w > 1e-5) ];
      wOp = wOp /sum(wOp);
	  wxOp = wOp || xOp;
   finish;

   use Grid; read all var {x1 x2} into xCan;

   alpha = 0.318; call AlphaOpt; create Aug318 var {w318 x1 x2}; append from wxOp;
   alpha = 0.531; call AlphaOpt; create Aug531 var {w531 x1 x2}; append from wxOp;
   alpha = 0.647; call AlphaOpt; create Aug647 var {w647 x1 x2}; append from wxOp;
   alpha = 0.949; call AlphaOpt; create Aug949 var {w949 x1 x2}; append from wxOp;


title2 "Optimal frequencies for 1000 runs";
proc sort data=f318; by x2 x1;
proc sort data=f531; by x2 x1;
proc sort data=f647; by x2 x1;
proc sort data=f949; by x2 x1;
data tab2001x; merge   f318   f531   f647   f949; by x2 x1;
   f318 = f318 / 100;
   f531 = f531 / 100;
   f647 = f647 / 100;
   f949 = f949 / 100;
   label f318 = "Alpha=0.318";
   label f531 = "Alpha=0.531";
   label f647 = "Alpha=0.647";
   label f949 = "Alpha=0.949";
proc print data=tab2001x noobs label;
   var x1 x2 f318 f531 f647 f949;
   format f318 f531 f647 f949 5.3;
run;
title2;

title2 "Optimal weights";
proc sort data=Aug318; by x2 x1;
proc sort data=Aug531; by x2 x1;
proc sort data=Aug647; by x2 x1;
proc sort data=Aug949; by x2 x1;
data tab2001c; merge Aug318 Aug531 Aug647 Aug949; by x2 x1;
   label w318 = "Alpha=0.318";
   label w531 = "Alpha=0.531";
   label w647 = "Alpha=0.647";
   label w949 = "Alpha=0.949";
proc print data=tab2001c noobs label;
   var x1 x2 w318 w531 w647 w949;
   format w318 w531 w647 w949 5.3;
run;
title2;

title2 "SAS Task 20.2 - Optimum 10-run design";
%Exact(10,0.318,f318,f318);
proc sort data=f318(rename=(COUNT=Reps)); by x2 x1;
data SASTask2002; merge f318 Aug318(rename=(w318=Wgt));
   RoundWgt = round(10*Wgt);
proc print data=SASTask2002 noobs;
   var x1 x2 Reps Wgt RoundWgt;
run;
title2;

title1;