/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 16.8                                              */
/*   TITLE: Fig. 16.9-16.11 - Optimally blocked and orthogonally      */
/*          blocked mixture designs                                   */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Uses the ADX macros to create simplex lattice designs.    */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Fig. 16.9-16.11 - Optimally blocked and orthogonally blocked mixture designs";


%adxgen;
%adxmix;

/*
/  The orthogonally blocked design of Fig. 16.9
/---------------------------------------------------------------------*/
data Orth;
   array x{3};
   input Block x1 x2 x3;
   do i = 1 to 3; x{i} = x{i} / 6; end; drop i;
   cards;
1 4 1 1
1 1 4 1
1 1 1 4
1 2 2 2
2 3 3 0
2 3 0 3
2 0 3 3
2 2 2 2
;

/*
/  The optimally blocked design from SLD(2), Fig. 16.10
/---------------------------------------------------------------------*/
%adxsld(Candidates,x1 x2 x3,2);
proc optex data=Candidates coding=orthcan seed=1;
   model x1|x2|x3@2 / noint;
   block structure=(2)4;
   output out=Des1;
run;

/*
/  The optimally blocked design from SLD(100), Fig. 16.11
/---------------------------------------------------------------------*/
%adxsld(Candidates,x1 x2 x3,10);
proc optex data=Candidates coding=orthcan seed=1;
   model x1|x2|x3@2 / noint;
   block structure=(2)4;
   output out=Des2;
run;

/*
/  To compute various criteria ...
/---------------------------------------------------------------------*/
proc iml;

/*
/  ... read in the designs, ...
/---------------------------------------------------------------------*/
   use Orth; read all var {x1 x2 x3} into fOrth;
             read all var {Block}    into bOrth;
   use Des1; read all var {x1 x2 x3} into fDes1;
             read all var {Block}    into bDes1;
   use Des2; read all var {x1 x2 x3} into fDes2;
             read all var {Block}    into bDes2;

/*
/  ... construct the second-order terms and ...
/---------------------------------------------------------------------*/
   do i = 1 to 2; do j = i+1 to 3;
      fOrth = fOrth || fOrth[,i]#fOrth[,j];
      fDes1 = fDes1 || fDes1[,i]#fDes1[,j];
      fDes2 = fDes2 || fDes2[,i]#fDes2[,j];
      end; end;

/*
/  ... the block terms, ...
/---------------------------------------------------------------------*/
   zOrth = designf(bOrth);
   zDes1 = designf(bDes1);
   zDes2 = designf(bDes2);
  
/*
/  ... combine them to form the entire model matrix, ...
/---------------------------------------------------------------------*/
   xOrth = zOrth || fOrth;
   xDes1 = zDes1 || fDes1;
   xDes2 = zDes2 || fDes2;

/*
/  ... print the orthogonality measure, and ...
/---------------------------------------------------------------------*/
   p = ncol(fDes1);
   print ((det(fOrth`*fOrth)*det(zOrth`*zOrth)/det(xOrth`*xOrth))**(-1/p))
         ((det(fDes1`*fDes1)*det(zDes1`*zDes1)/det(xDes1`*xDes1))**(-1/p))
         ((det(fDes2`*fDes2)*det(zDes2`*zDes2)/det(xDes2`*xDes2))**(-1/p));

/*
/  ... compute D and relative D criteria.
/---------------------------------------------------------------------*/
   dOrth = det(xOrth`*xOrth)**(1/ncol(xOrth));
   dDes1 = det(xDes1`*xDes1)**(1/ncol(xDes1));
   dDes2 = det(xDes2`*xDes2)**(1/ncol(xDes2));
   print dOrth  dDes1         dDes2       ,
         .     (dDes1/dOrth) (dDes2/dDes1);
quit;
title1;