/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 14.2                                              */
/*   TITLE: Table 14.1 - D-optimum exact N-trial designs for a        */
/*          second-order model in m quantitative factors with one     */
/*          qualitative factor at B levels with one qualitative       */
/*          factor                                                    */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Took ~3 minutes to run on server-class PC in May, 2007.   */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Table 14.1 - D-optimum exact second-order designs";
title2 "with quantitative factors";

/*
/  Table 14.1 - A version of AD Table 13.1, but without fixed block
/  sizes.
/---------------------------------------------------------------------*/
%macro iBQuad(m,B,n,nRandom=10);

options nonotes;

/*
/  Create a macro variable that holds the quadratic model.
/---------------------------------------------------------------------*/
%let xmodel = x1;                /* Linear effects and cross-products */
%do j = 2 %to &m; %let xmodel = &xmodel|x&j; %end;
%let xmodel = &xmodel@2;
                                 /* quadratic effects                 */
%do j = 1 %to &m; %let xmodel = &xmodel x&j*x&j; %end;

/*
/  Create candidate points (-1,0,1)**m and add "standard order" IDs.
/---------------------------------------------------------------------*/
proc factex;
   factors x1-x&m / nlev=3;
   output out=Quantitative;
proc sort data=Quantitative;
   by x&m-x1;
data Quantitative; set Quantitative;
   id = _N_;
proc factex;
   factors Block / nlev=&B;
   output out=Candidates
      Block nvals=(1 to &B)
      pointrep=Quantitative;
run;

/*
/  Find the D-optimal design with OPTEX and save the efficiencies.
/---------------------------------------------------------------------*/
proc optex data=Candidates coding=orthcan noprint;
   class Block;
   model Block &xmodel;
   generate n=&N method=m_fedorov niter=100000 keep=1 initdesign=random;
   id id;
   output out=opDesign;
run;

%let n1 = .;
%let n2 = .;
%let n3 = .;
proc freq data=opDesign noprint;
   table Block/ out=BSize;
proc sort data=BSize;
   by Count;
data _null_; set BSize;
   call symput('n'||trim(left(_N_)),trim(left(Count)));
run;

proc factex;
   factors x1-x&m / nlev=5;
   output out=Gridx
      %do i = 1 %to &m;
         x&i nvals=(-1 to 1 by .5)
         %end;
      ;
run;
   factors Block / nlev=&B;
   output out=Grid
      Block nvals=(1 to &B)
      pointrep=Gridx;
run;

/*
/  Now, we use non-linear optimization in IML both to improve the
/  exact design and to find the optimal design weights.  We begin by
/  defining functions that will be required.
/---------------------------------------------------------------------*/
proc iml;

   /*
   /  From given blocks B and quantitative factors x, create the
   /  extended design matrix.
   /------------------------------------------------------------------*/
   start MakeZ(B,x);
      Z = design(B);                  /* Indicators of levels of B    */
      do i = 1 to ncol(x);            /* Linear terms in x            */
         Z = Z || x[,i];
         end;
      do i = 1 to ncol(x)-1;          /* Cross-products in x          */
         do j = i+1 to ncol(x);
            Z = Z || x[,i]#x[,j];
            end;
         end;
      do i = 1 to ncol(x);            /* Quadratic terms in x         */
         Z = Z || x[,i]#x[,i];
         end;
      return(Z);
   finish;

   /*
   /  From given values xx of the quantitative factors and for fixed
   /  values B of the blocks, compute the exact design's D-criterion
   /  log(|Z'Z|**(1/p)).
   /------------------------------------------------------------------*/
   start DCrit(xx) global(B);
      x = shape(xx,nrow(xx)*ncol(xx)/&m,&m);
      Z = MakeZ(B,x);
      d = det(Z`*Z/nrow(Z));
      if (d > 0) then return(log(d)/ncol(Z));
      else            return(-9999);
   finish;

   /*
   /  From given weights w and for a fixed candidate design matrix Z,
   /  compute the weighted design's D-criterion log(|Z'WZ|**(1/p)).
   /------------------------------------------------------------------*/
   start DCritW(w) global(Z);
      Dw = diag((w##2)/sum(w##2));   /* Use relative squared weights,
                                        so we don't need constraints. */
      return(log(det(Z`*Dw*Z)**(1/ncol(Z))));
   finish;

/*
/  Improve the OPTEX-generated design using constrained non-linear
/  optimization to adjust the quantitative factor values so that the
/  the exact design's D-criterion is optimized.
/---------------------------------------------------------------------*/
   use opDesign;
   read all var ("x1":"x&m") into xop;
   read all var ("Block")    into B;

   xinit = shape(xop,&m*nrow(xop),1); /* Initialize optimization at
                                         OPTEX's design               */
   con =    shape(-1,1,nrow(xinit))   /* Lower and upper constraints  */
         // shape( 1,1,nrow(xinit));
   call nlpqn(rc,xx,"DCrit",xinit,{1},con);
   xim = shape(xx,nrow(xx)*ncol(xx)/&m,&m);

/*
/  It's also easy to try optimizing DCrit directly, from random
/  initial designs.
/---------------------------------------------------------------------*/
   dmax = DCrit(xx);
   do i = 1 to &nRandom;
      x1 = shape(2*ranuni(1:nrow(xinit)) - 1,nrow(xinit),1);
      call nlpqn(rc,xx,"DCrit",x1,{1},con);
      if (DCrit(xx) > dmax) then do;
         xim = shape(xx,nrow(xx)*ncol(xx)/&m,&m);
         dmax = DCrit(xx);
         end;
      end;
   
/*
/  Save the optimum values of the D-criterion.  Take the p-th power
/  since that's what AD table.
/---------------------------------------------------------------------*/
   p = &B + &m + (&m*(&m-1)/2) + &m;
   call symput("DOp",trim(left(char(exp(p*DCrit(xop)),20))));
   call symput("DIm",trim(left(char(exp(p*DCrit(xim)),20))));

/*
/  For D-efficiency calculations, we need the optimal design weights.
/  This is also a job for non-linear optimization, but now on the
/  design weights over a fixed candidate Grid.  Optimize the weights
/  and save the maximal weighted D-criterion.
/---------------------------------------------------------------------*/
   use Grid;
   read all var ("x1":"x&m") into xCan;
   read all var ("Block")    into BCan;
   Z = MakeZ(BCan,xCan);

   winit = j(nrow(xCan),1) / nrow(xCan);   /* Initialize with uniform */
   call nlpqn(rc,w,"DCritW",winit,1);
   call symput('DOpt',trim(left(char(exp(DCritW(w)),20))));

/*
/  Accumulate design information in the data set Table14_1.
/---------------------------------------------------------------------*/
data Eff;
   m = &m;
   B = &B;
   n1 = &n1;
   n2 = &n2;
   n3 = &n3;
   p  = B + m + (m*(m-1)/2) + m;
   N = &N;
   DOp = &DOp;
   DIm = &DIm;
   DOpt = &DOpt;
data Table14_1; set Table14_1 Eff;
run;
options notes;

%mend;

data Table14_1; if (0); run;
%iBQuad(1,2,4 );
%iBQuad(1,2,5 );
%iBQuad(1,2,6 );
%iBQuad(1,3,5 );
%iBQuad(1,3,6 );
%iBQuad(1,3,7 );
%iBQuad(1,3,8 );
%iBQuad(1,3,9 );
              
%iBQuad(2,2,7 );
%iBQuad(2,2,8 );
%iBQuad(2,2,9 );
%iBQuad(2,2,10);
%iBQuad(2,2,11);
%iBQuad(2,2,12);
%iBQuad(2,2,13);
%iBQuad(2,2,14);
%iBQuad(2,2,15);
%iBQuad(2,2,16);
%iBQuad(2,2,17);
%iBQuad(2,2,18);
              
%iBQuad(2,3,8 );
%iBQuad(2,3,9 );
%iBQuad(2,3,10);
%iBQuad(2,3,11);
%iBQuad(2,3,12);
%iBQuad(2,3,15);
%iBQuad(2,3,18);

/*
/  The efficiencies are functions of the the p-th roots of the D-
/  criteria; the criterion DOpt for the optimal design weight is
/  already the p-th root.
/---------------------------------------------------------------------*/
data Table14_1; set Table14_1;
   DEff   = ((DOp**(1/p))) / DOpt;
   DEffAA = ((DIm**(1/p))) / DOpt;
proc print data=Table14_1(drop=DOpt) noobs;
   var m B N p n1 n2 n3 DOp DIm DEff DEffAA;
   format DOp DIm     e9.3;
   format DEff DEffAA 6.4;
run;

title2;
title1;