/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 18.2                                              */
/*   TITLE: Example 18.2: Exponential decay                           */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 18.2: Exponential decay";

title2 "Table 18.4 - Comparisons of optimum designs";
proc iml;

   /*
   /  This function simplifies a design, often revealing an equivalent
   /  design with fewer support and thus one whose structure is easier
   /  to see.
   /------------------------------------------------------------------*/
   start Consolidate(x,w,ex,ew);

      x = x[loc(w > ew)]`;             /* Drop points with small wgts */
      w = w[loc(w > ew)]`;

      Group = shape(0,1,ncol(x));      /* Replace nearby points with  */
      do i = 1 to ncol(x);             /* their avg, with wgts equal  */
         if (^Group[i]) then do;       /* to sum(wgt) over averaged   */
            Group[i] = max(Group) + 1; /* points.                     */
            do j = i+1 to ncol(x);
               if (abs(x[i] - x[j]) < ex) then
                  Group[j] = Group[i];
               end;
            end;
         end;
      xNew = shape(.,1,max(Group));
      wNew = shape(.,1,max(Group));
      do i = 1 to max(Group);
         xNew[i] = (x[loc(Group = i)])[:];
         wNew[i] = (w[loc(Group = i)])[+];
         end;

      r = rank(xNew);                  /* Make sure the points are    */
      x = xNew;                        /* sorted by x.                */
      w = wNew;
      x[r] = xNew;
      w[r] = wNew;
   finish;

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

   tPop = 1/7 || 1/sqrt(7) || 1 || sqrt(7) || 7;
   pPop = j(1,ncol(tPop))/ncol(tPop);

   start DCrit(w,x) global(tPop,pPop,Criterion);
      if      (Criterion = 1) then EldMi = 0;
      else if (Criterion = 2) then E_dMi = 0;
      else if (Criterion = 3) then E__Mi = shape(0,1,1);
      else if (Criterion = 4) then E_dM_ = 0;
      else                         E__M_ = shape(0,1,1);
      do i = 1 to ncol(tPop);
         fi = x#exp(-tPop[i]*x);
         Mi = fi*diag(w)*fi`;
         if      (Criterion = 1) then EldMi = EldMi + pPop[i]*llog(det(ginv(Mi)));
         else if (Criterion = 2) then E_dMi = E_dMi + pPop[i]*     det(ginv(Mi));
         else if (Criterion = 3) then E__Mi = E__Mi + pPop[i]*         ginv(Mi) ;
         else if (Criterion = 4) then E_dM_ = E_dM_ + pPop[i]*     det(     Mi );
         else                         E__M_ = E__M_ + pPop[i]*              Mi  ;
         end;
      if      (Criterion = 1) then return(-         EldMi  );
      else if (Criterion = 2) then return(-llog(    E_dMi ));
      else if (Criterion = 3) then return(-llog(det(E__Mi)));
      else if (Criterion = 4) then return( llog(    E_dM_ ));
      else                         return( llog(det(E__M_)));
   finish;

   start DCritW(w) global(x);
      return(DCrit((w##2)/sum(w##2),x));
   finish;

   start DCritX(x) global(w);
      return(DCrit(w,x));
   finish;

   start DCritWX(wx);
      w = shape(wx,2,ncol(wx)/2)[1,];
      w = (w##2)/sum(w##2);
      x = shape(wx,2,ncol(wx)/2)[2,];
      return(DCrit(w,x));
   finish;

   /*
   / 
   /------------------------------------------------------------------*/
   options nonotes;
   do Criterion = 1 to 5;

      nint = 50;
      x    = 10*(1:nint)/nint;
      winit = j(ncol(x),1) / ncol(x);
      call nlpqn(rc,wGrid,"DCritW",winit,1);
      wGrid = ((wGrid##2)/sum(wGrid##2));
      xGrid = x  [loc(wGrid > 0)];
      call Consolidate(xGrid,wGrid,0,1e-4);
      xOp = xGrid;
      wOp = wGrid;

      i = 0;
      Improvement = 1;
      do while ((i < 100) & (Improvement > 1e-5));
         i = i + 1;

         last = (xOp` || wOp`);

         wxinit = sqrt(wOp) || xOp;
         con =    (shape(.,1,ncol(xOp)) || shape( 0,1,ncol(xOp)))
               // (shape(.,1,ncol(xOp)) || shape(10,1,ncol(xOp)));
         call nlpqn(rc,wxOp,"DCritWX",wxinit,1,con);
         wOp = shape(wxOp,2,ncol(wxOp)/2)[1,];
         xOp = shape(wxOp,2,ncol(wxOp)/2)[2,];
         wOp = (wOp##2)/sum(wOp##2);
         call Consolidate(xOp,wOp,1e-2,1e-4);

         if (ncol(xOp) ^= nrow(last)) then
            Improvement = 999;
         else
            Improvement = max(abs(last - (xOp` || wOp`)));
         end;

      Support = xOp`;
      Weights = wOp`;
      print Criterion Support [format=6.4] Weights [format=6.4];
      end;

title2;
title1;