/*
   SAS/IML program for calculating critical values of multiple comparison procedures.

   Author : Frank Bretz and Paul Somerville
   Contact: bretz@ifgb.uni-hannover.de
   Date   : 11.05.99 - 1.0 - start version
            20.02.00 - 1.1 - inclusion of several new swtypes (7-15 and 20) and some minor changes
            02.06.00 - 1.2 - inclusion of several new swtypes (7-15); based on version 1.0 (modifications of version 1.1 have been dropped)

   Input  : (Required parameters)
            swcov   : Flag to include var-covar matrix (=0) or not (=1)
            swtype  : Type of multiple comparison procedure
            conf    : Confidence level
            givense : User requested SE of estimates, if alternate procedure is used
            k       : Number of treatment groups
            seed    : Seed for random number generator
            ndenom  : Degrees of freedom
            mocar   : Number of random directions
            irep    : Number of prliminary estimates (if negative, alternate procedure is performed)

            (Optional parameters)
            mm      : Number of contrasts (if swtype = 0)
            vc      : Var-covar matrix (if swcov = 0)
            contr   : Contrasts (if swtype = 0 or 20)
            icontrol: Control group {if [(swtype=4 OR 5) AND (swcov=0)] OR [(swtype=9)]}
            iflo    : Lower bound for peak (swtype = 10)
            ighi    : Higher bound for peak (swtype = 10)
            znn     : Vector of sample sizes (swtype = 11-15)

   Output :

   References: Somerville (1997, 1998, 1999), Somerville and Bretz (2000, 2001)
*/

options nosource;

proc iml;


  /* INSERT SITUATION PARAMETERS HERE */


  swcov  = 0     ;
  swtype = 15     ;
  conf   = .975    ;
  givense= 999   ;
  k      = 4     ;
  seed   = 1893  ;
  ndenom = -1    ;
  mocar  = 1000  ;
  irep   = 10    ;
  mm     = 99    ;
  vc={110.25     .    .     .,
      130.2    255.77    .    .,
      139.65    318.44    526.74   .,
      149.1    340.71    630.64    898.66};
  contr=.;
  icontrol=.;
  znn={10 11 12 13};



  /* END SITUATION PARAMETERS */






  if ncol(icontrol)=0 then icontrol=.;
  if ncol(vc)=0 then vc=.;
  if ncol(contr)=0 then contr=.;
  nq     = 25    ;
  if swcov=0 then do;
    do i=1 to k-1;
      do j=i+1 to k;
        vc[i,j]=vc[j,i];
      end;
    end;
  end;

/***********************************************************
 ***********************************************************
 ***********************************************************/

      /* This function gives the logarithm of the gamma function when
         the argument is a positive integer or an integer divided by 2 */

      start gln(a);
        i1=int(a);
        c2=0;
        c1=1;
        gln=0;
        if (a-i1)=0 then i1=i1-2;
        else do;
          c1=-.5;
          gln=log(sqrt(arcos(-1)));
        end;
        do  i=1 to i1;
         z=c1+i;
         x=log(z);
         gln=gln+x;
        end;
        return(gln);
      finish;


      /* This function gives the ratio of gamma(a+b) to gamma(a)*gamma(b) */

      start beta(a,b);
        beta=gln(a+b)-gln(a)-gln(b);
        beta= exp(beta);
        return(beta);
      finish;


      /* Calculates k-1 rows of Helmert matrix */

      start helmert(k,k1);
        hm=j(k1,k,.);
        do i=1 to k1;
          zi=i;
          do j=k-i to k;
            hm[i,j]=1/sqrt(zi*(zi+1));
          end;
          do j=1 to k1-i;
            hm[i,j]=0;
          end;
          hm[i,k-i]=-zi*hm[i,k-i];
        end;
        return(hm);
      finish;


      /* SUBROUTINE TO ENUMERATE ALL OF THE PAIRS OF FACES */

      start facelist(mm) global(k, ifaceb,ifacee);
        ifacee=j(1,mm,.);
        ifaceb=j(1,mm,.);
        i=0;
        ie=2;

        do until(ie=k+1);
          ib=1;
          krit=0;
          do until(krit=1);
            i=i+1;
            ifaceb[i]=ib;
            ifacee[i]=ie;
            if (ib^=ie-1) then ib=ib+1;
            else krit=1;
          end;
          ie=ie+1;
        end;
      finish;


      /* Generation of Contrast Matrices */

      start contrast(swtype,mm,icontrol) global(znn,ighi,iflo,k,k1,ifaceb,ifacee);
        contr=j(mm,k,0);
        if (swtype=1 | swtype=2 | swtype=6) then run facelist(mm);

        if swtype=1 then do;
          imm=mm/2;
          do i=1 to imm;
            j=i+imm;
            contr[i,ifaceb[i]]=1;
            contr[i,ifacee[i]]=-1;
            contr[j,ifaceb[i]]=-1;
            contr[j,ifacee[i]]=1;
          end;
        end;

        if (swtype=2 | swtype=6) then do;
          do i=1 to mm;
            contr[i,ifaceb[i]]=1;
            contr[i,ifacee[i]]=-1;
          end;
        end;

        if (swtype=3 | swtype=4) then do;
          do i=1 to mm;
            contr[i,icontrol]=-1;
          end;

          do i=1 to icontrol-1;
            contr[i,i]=1;
          end;

          do i=icontrol+1 to mm;
            contr[i,i+1]=1;
          end;

          if icontrol<=mm then contr[icontrol,icontrol+1]=1;
        end;

        if swtype=5 then do;
          do i=1 to k1;
            j=i+k1;
            contr[i,icontrol]=-1;
            contr[j,icontrol]=1;
          end;

          do i=1 to icontrol-1;
            j=i+k1;
            contr[i,i]=1;
            contr[j,i]=-1;
          end;

          do i=icontrol+1 to k1;
            j=i+k1;
            contr[i,i+1]=1;
            contr[j,i+1]=-1;
          end;

          if icontrol<=k1 then do;
            contr[icontrol,icontrol+1]=1;
            contr[icontrol+k1,icontrol+1]=-1;
          end;
        end;

        if (swtype=7 | swtype=8) then do;
          do i=1 to k-1;
            contr[i,i]=1;
            contr[i,i+1]=-1;
          end;
        end;

        if (swtype=8) then do;
          do i=k to mm;
            contr[i,i-k+1]=-1;
            contr[i,i-k+2]=1;
          end;
        end;

        if (swtype=9) then do;
          mm1=icontrol*(icontrol-1)/2;
          mm2=(k+1-icontrol)*(k-icontrol)/2;
          if (mm1 ^= 0) then do;
            ifacee=j(1,mm,0);
            ifaceb=j(1,mm,0);
            i=0;
            ie=2;
            do until(ie=icontrol+1);
              ib=1;
              do until(ib=ie);
                i=i+1;
                ifaceb[i]=ib;
                ifacee[i]=ie;
                ib=ib+1;
              end;
              ie=ie+1;
            end;
            do jj=1 to mm1;
              contr[jj,ifaceb[jj]]=-1;
              contr[jj,ifacee[jj]]=1;
            end;
          end;

          if (mm2 ^= 0) then do;
            ifacee=j(1,mm,0);
            ifaceb=j(1,mm,0);
            i=0;
            ie=2;
            do until(ie=k-icontrol+2);
              ib=1;
              do until(ib=ie);
                i=i+1;
                ifaceb[i]=ib;
                ifacee[i]=ie;
                ib=ib+1;
              end;
              ie=ie+1;
            end;
            do jj=1 to mm2;
              contr[jj+mm1,ifaceb[jj]+icontrol-1]=1;
              contr[jj+mm1,ifacee[jj]+icontrol-1]=-1;
            end;
          end;
        end;

        if (swtype=10) then do;
          in=0;
          do i=1 to ighi-1;
            do j=i+1 to ighi;
              in=in+1;
              contr[in,i]=1;
              contr[in,j]=-1;
            end;
          end;
          do j=iflo to k-1;
            do i=j+1 to k;
              in=in+1;
              contr[in,j]=-1;
              contr[in,i]=1;
            end;
          end;
        end;

        if swtype=11 then contr=cm_will(znn);
        if swtype=12 then contr=cm_marc(znn);
        if swtype=13 then contr=cm_hirot(znn);
        if swtype=14 then contr=cm_mcder(znn);
        if swtype=15 then contr=cm_iso(znn);


print contr;
        return(contr);
      finish;


  start cm_will(n);
    k1=ncol(n)-1;

    c=j(1,k1,.);
    cm=j(k1,k1,0);
    cm=j(k1,1,-1)||cm;

    do i=1 to k1;
      x=sum(n[k1-i+2:k1+1]);
      do j=1 to i;
        cm[i,k1-j+2]=n[k1-j+2]/x;
      end;
    end;

    return(cm);
  finish;

  start cm_marc(n);
    k=ncol(n);

    cm1=j(k-1,k,0);
    cm2=cm1;

    do i=1 to k-1;
      cm1[i,i+1:k]=t(n[i+1:k]/sum(n[i+1:k]));
    end;

    do i=1 to k-1;
      cm2[i,1:i]=t(n[1:i]/sum(n[1:i]));
    end;

    row=k*(k-1)/2;
    cm=j(row,k,0);
    index=1;
    do i=1 to k-1;
      do j=1 to i;
        cm[index,]=cm1[i,]-cm2[j,];
        index=index+1;
      end;
    end;

    return(cm);
  finish;

  start cm_hirot(n);
    k=ncol(n);

    cm=j(k-1,k,0);
    do i=1 to k-1;
      do j=1 to k;
        if j>i then cm[i,j]=n[j]/sum(n[i+1:k]);
        else cm[i,j]=-n[j]/sum(n[1:i]);
      end;
    end;

    return(cm);
  finish;

  start cm_mcder(n);
    k=ncol(n);

    cm=j(k-1,k,0);
    do i=1 to k-1;
      do j=1 to k;
        if j <= i then cm[i,j]=-n[j]/sum(n[1:i]);
        else if j=i+1 then cm[i,j]=1;
      end;
    end;

    return(cm);
  finish;


  /******************************************************************/
  /* The module SCHNITT cuts the place STELLE out of the vector VEC */
  /* and returns a vector which dimension dimension is reduced by 1 */
  /******************************************************************/


  start schnitt(vec,stelle);
    col=ncol(vec);
    if stelle=col then neu=vec[,1:col-1];
    else do;
      neu=vec[,1:col-1];
      neu[stelle:col-1]=vec[,stelle+1:col];
    end;
    return(neu);
  finish;


  /******************************************************************/
  /* The main module to calculate the contrast coefficients         */
  /******************************************************************/


  start cm_iso(n);
    k=ncol(n)-1;
    anz=2##k;
    cm=j(2##k-1,k+1,.);
    mat=j(anz,k,0);
    mat[2,k]=1;


    /*****************************************/
    /* Generation of all 2^k-1 possibilities */
    /*   0: two subsequent groups are pooled */
    /*   1: subsequent groups are not pooled */
    /*****************************************/


    do i=2 to k;
      mat[2##(i-1)+1:2##i,k-i+1]=j(2##(i-1),1,1);
      mat[2##(i-1)+1:2##i,k-i+2:k]=mat[1:2##(i-1),k-i+2:k];
    end;
    mat=mat[2:2##k,];

    row=nrow(mat);

    n1=j(1,k,n[1]);
    n2=j(1,k,sum(n[2:k+1]));
    do j=2 to k;
      n1[j]=sum(n[1:j]);
      n2[j]=sum(n[j+1:k+1]);
    end;


    /*************************************************/
    /* Beginning of main loop: each row = one of the */
    /* 2^k-1 possibilities is evaluated on its own   */
    /*************************************************/


    do i=1 to row;
      y=mat[i,];                                    /* Amalgamated contrast vector */
      x=j(1,k,1);                                   /* Weight vector according to the number of amalgamations */

      s=mat[i,+];
      ssd=j(1,s+1,0);
      nn=-y#sqrt(n1#n2/n[+]);
      nn_help=nn;
      do j=1 to ncol(nn);
        nn[j]=nn_help[ncol(nn)-j+1];
      end;

      krit=0;
      j=1;
      do until(krit=1);
        if nn[j]=0 then nn=schnitt(nn,j);
        else j=j+1;
        if j>ncol(nn) then krit=1;
      end;

      ssd[2:s+1]=nn;
      ssd_help=ssd;
      do j=1 to ncol(ssd);
        ssd[j]=ssd_help[ncol(ssd)-j+1];
      end;


      krit=0;
      j=1;
      y_help=y;
      n_gew=n;

      do until(krit=1);
        if y_help[j]=0 then do;
          n_gew[j]=n_gew[j]+n_gew[j+1];
          n_gew=schnitt(n_gew,j+1);
          y_help=schnitt(y_help,j);
        end;
        else j=j+1;
        if all(y_help)=1 then krit=1;
      end;


      krit=0;
      j=1;

      do until(krit=1);                             /* Actual amalgamation */
        if y[j]=0 then do;
          if y[j]=y[j+1] then do;
            y=schnitt(y,j+1);                       /* Pooling, i.e. reduction of contrast vector */
            x=schnitt(x,j+1);                       /* Reduction des weight vector */
            x[j]=x[j]+1;                            /* Increase of weight by 1, if pooling happend */
          end;
          else j=j+1;
        end;
        else j=j+1;
        if j>=ncol(y) then krit=1;
      end;


      do j=1 to ncol(y);                            /* Correction of weights */
        if y[j]=0 then x[j]=x[j]+1;
      end;

      j=1;
      do while(j<(ncol(y)-1));                      /* Special case. if a non-pooling group lies between */
        if y[j]=0 then do;                          /*    two groups to be pooled, e.g. 010              */
          l=j+1;
          do until(y[l]=0 | l>=ncol(y));
            l=l+1;
          end;
          if y[l]=0 then do;
            y=schnitt(y,j+1);
            x=schnitt(x,j+1);
          end;
        end;
        j=j+1;
      end;

      if i=row then do;                             /* Special case, if no pooling at all   */
        y=y||{1};                                   /* (always the last of the 2^k-1 cases) */
        x=x||{1};
      end;

      c=j(1,s+1,ssd[1]/n_gew[1]);
      do j=2 to s+1;
        c[j]=(ssd[j]-sum(c[1:j-1]#n_gew[1:j-1]))/n_gew[j];
      end;

      contrast=c;
      cm[i,1:x[1]]=contrast[1];                     /* Solution of amalgamation */
      do j=2 to ncol(y);
        cm[i,sum(x[1:j-1])+1:sum(x[1:j])]=contrast[j];
      end;
    end;
    cm=cm#n;

    return(cm);
  finish;





/***********************************************************
 ***********************************************************
 ***********************************************************/


      /* This routine does the binning of the reciprocal distances */

     start fillbin2(recdist,nr) global(pl, ibin);
       i=nr;
       krit=0;
       do until(krit=1 | i=0);
         if recdist>= pl[i] then do;
           j=nr+1-i;
           ibin[j]=ibin[j]+1;
           krit=1;
         end;
       else i=i-1;
       end;
     finish;


     /* This routine calculates proportions in the bins */

     start binprop(points,nq) global(ibin, prop);
       do  i=2 to nq;
         ibin[i]=ibin[i-1]+ibin[i];
       end;
       do  i=1 to nq;
         prop[i]=1-ibin[i]/points;
       end;
     finish;


     /* This routine calculates the probability of correct statements outside of region of radius q. */

     start fringe(nr,prop) global(wl,pl,zm,zn,fconst);
       x=0;
       do i=1 to nr;
         j=nr-i+1;
         cc=zm*log(1+zn*pl[i]*pl[i])/2;
         cd=zn*log(1+1/zn/pl[i]/pl[i])/2;
         xx=exp(-log(pl[i])-cc-cd);
         x=x+xx*wl[i]*prop[j];
       end;
       fringe=x*fconst;  *print wl pl fconst;
       return(fringe);
     finish;


     /* Function obtains confidence corresponding to values of q, k1 and ndenom */

     start prob(q,nq,prop) global(zm,zn);

       qq=sqrt(2)/q;
       run gauleg(0,qq,nq);
       if(zn>0) then do;
         x=fringe(nq,prop);
         qqdk1=q*q/zm/2;
         poff=probf(qqdk1,zm,zn);
       end;
       else do;
         x=fringinf(nq,zm,prop);
         qqdk1=q*q/2;
         poff=probchi(qqdk1,zm);
       end;
       prob2=x+poff;
       return(prob2);
     finish;


     /* Calculation of g(v) for Multivariate Normal Distribution */

     start fringinf(nr,q,prop) global(wl,pl);
       x=0;
       do i=1 to nr;
         j=nr-i+1;
         ce=-log(gamma(q/2))-(q/2-1)*log(2)-(q+1)*log(pl[i])-1/pl[i]/pl[i]/2;
         xx=exp(ce);
         x=x+xx*wl[i]*prop[j];
       end;
       fringinf=x;
       return(fringinf);
     finish;


     /* Determination of Nodes and Weigths for Gauss-Legendre */

     start gauleg(x1,x2,n) global(pl,wl);
       x=j(1,n,.);
       w=j(1,n,.);

       eps2=1E-09;
       m=(n+1)/2;
       xm=(x2+x1)/2;
       xl=(x2-x1)/2;
       pi=2*arsin(1);
       do i=1 to m;
         z=cos(pi*(i-0.25)/(n+0.5));
         do until (abs(z-z1)<eps2);
           p1=1;
           p2=0;
           do j=1 to n;
             p3=p2;
             p2=p1;
             p1=((2*j-1)*z*p2-(j-1)*p3)/j;
           end;
           pp=n*(z*p1-p2)/(z*z-1);
           z1=z;
           z=z1-p1/pp;
         end;
         x[i]=xm-xl*z;
         x[n+1-i]=xm+xl*z;
         w[i]=2*xl/((1-z*z)*pp*pp);
         w[n+1-i]=w[i];
       end;

       pl=x;
       wl=w;
     finish;


     /* INES according to Ridders */

     start ridders(prop,nq) global(brakl,braku,conf,xacc);
       x1=brakl;
       x2=braku;
       alpha=conf;
       quan_old=8;
       index=0;

       f1=prob(x1,nq,prop)-alpha;
       f2=prob(x2,nq,prop)-alpha;

       do until(abs(diff)<xacc);
         x=(x1+x2)/2;
         f=prob(x,nq,prop)-alpha;

         wurz=sqrt(f*f-f1*f2);
         quan=x+(x-x1)*sign(f1-f2)*f/wurz;
         fquan=prob(quan,nq,prop)-alpha;
         diff=quan-quan_old;
         quan_old=quan;

         if f*fquan<0 then do;
           x1=x;
           f1=f;
           x2=quan;
           f2=fquan;
         end;
         else do;
           if f1*fquan<0 then do;
             x2=quan;
             f2=fquan;
           end;
           else do;
             if f2*fquan<0 then do;
               x1=quan;
               f1=fquan;
             end;
           end;
         end;

         index=index+1;
       end;
       return(quan);
     finish;


     /* This subroutine calculate the distances from the origin to the boundaries */

     start mvidist(dc,gg,k1) global(mm, dis1,dis2);
       dis1=0;
       dis2=0;
       do i=1 to mm;
         a=0;
         do j=1 to k1;
          a=a+dc[j]*gg[i,j];
         end;
         if a<0 then dis1=min(a,dis1);
         else dis2=max(a,dis2);
       end;
       dis1=-dis1;
     finish;



     start qvalue(nq,mocar) global(ibin,pl,wl,dis1,dis2,k1,irr,gg,prop);

       run gauleg(0,1,nq);

       ibin=j(1,nq,0);
       do index2=1 to mocar;
         dc=rannor(j(k1,1,irr));
         dc=dc/sqrt(sum(dc##2));

         run mvidist(dc,gg,k1);
         run fillbin2(dis1,nq);
         run fillbin2(dis2,nq);
       end;

       points=mocar*2;
       run binprop(points,nq);
*print points prop;
       qqq=ridders(prop,nq);
       return(qqq);
     finish;




      start estimate(i,nq,j) global(k,givense,sx,zrep,zmocar,irep,qans, mocar,zmean,se);
        beta={-6.55576 -8.1885 -8.2998 -7.5002 -7.5002 -7.2686 -8.2998,
                .20617   .19807  .19615  .19964  .19964  .16873  .19615,
              -3.8513  -3.9409 -3.7847 -2.9872 -2.9872 -3.4967 -3.7847};
        beta=t(beta);

        if j=-1 then do;
          mocar = 10000;
          coc=qvalue(nq,mocar);
          zk=k;
          newmocar=1000000/(givense)**2*(exp(beta[i+1,1]+beta[i+1,2]*coc+beta[i+1,2]/zk))**2;
          add = newmocar-mocar;
          if add>0 then do;
            tt=coc;
            mocar=int(max(add,256));
            coc=qvalue(nq,mocar);
            coc=(10000*tt+mocar*coc)/(10000+mocar);
            mocar=mocar+10000;
          end;
          zmean = coc;
          se=givense;
          if add<0 then se=exp(beta[i+1,1]+beta[i+1,2]*coc+beta[i+1,2]/zk)*10;
        end;

        if j=0 then do;
          coc=qvalue(nq,mocar);
          zmean=sx/zrep + coc;
          zk=k;
          se=exp(beta[i+1,1]+beta[i+1,2]*coc+beta[i+1,2]/zk)*1000/sqrt(zmocar);
        end;

        if j=1 then do;
          mocar=zmocar;
          sx=0;
          ssx=0;
          do njj=1 to irep;
            qqq=qvalue(nq,mocar);
*print qqq;
            qans[njj]=qqq;
            coc=qans[1];
            cm=qans[njj]-coc;
            sx=sx+cm;
            ssx=ssx+cm*cm;
          end;
          sd=(ssx-sx*sx/zrep)/(zrep-1);
          se=sqrt(sd/zrep);
          sd=sqrt(sd);
          zmean=sx/zrep+coc;
          mocar=zrep*mocar;
        end;
*        print zmean;
      finish;



/********************************************************************
 ********************************************************************
 ********************************************************************
 ********************************************************************/


      zmocar=mocar;
      if (swtype^=4 & swtype^=5 & swtype^=9)  then icontrol=1;
      if (swtype=1)              then mm=k*(k-1);
      if (swtype=2 | swtype=6)   then mm=k*(k-1)/2;
      if (swtype=3 | swtype=4)   then mm=k-1;
      if (swtype=5)              then mm=2*(k-1);
      if (swtype=7)                           then mm=k-1;
      if (swtype=8)                           then mm=2*(k-1);
      if (swtype=9)                           then do;
                                                     mm1=icontrol*(icontrol-1)/2;
                                                     mm2=(k+1-icontrol)*(k-icontrol)/2;
                                                     mm=mm1+mm2;
                                                   end;
      if (swtype=10)                          then mm=ighi*(ighi-1)/2+(k-iflo+1)*(k-iflo)/2;
      if (swtype=11)                          then mm=k-1;
      if (swtype=12)                          then mm=k*(k-1)/2;
      if (swtype=13)                          then mm=k-1;
      if (swtype=14)                          then mm=k-1;
      if (swtype=15)                          then mm=2**(k-1)-1;
      if (swtype^=0)                          then contr=j(mm,k,.);
      if (swtype^=10)                         then do;
                                                     iflo=.;
                                                     ighi=.;
                                                   end;



      brakl = .0001;
      braku = 5;
      xacc  = .0001;
      zero=0;
      zone=1;
      two=2;
      iseed=abs(seed);
      jseed=iseed;


      /* Random Number Generation */

      iy=0;
      iran=j(1,97,.);
      mr=714025;
      ia=1366;
      ic=150889;
      rm=1/mr;
      do  j=1 to 97;
        iseed=mod(ia*iseed+ic,mr);
        iran[j]=iseed;
      end;


      irr=iseed;
      k1=k-1;
      zm=k1;
      zn=ndenom;
      zmd2=zm/2;
      zmnd2=(zm+zn)/2;
      zmdzn=zm/zn;
      znd2=zn/two;
      if ndenom>0 then cnlogn=zn*log(zn)/2;
      sx=0;
      ssx=0;
      zrep=irep;



      dc=j(1,k1,.);

      iii=10*((irep-1)/10)+10;
      qans=j(1,iii,0);
      lines=(irep-1)/10+1;

      if ndenom>0 then fconst=2*beta(zmd2,znd2);
      if swcov^=0 then vc=I(k);
      hm=helmert(k,k1);
      vcc=hm*vc*t(hm);
      c=t(root(vcc))+1E-12;
      p=t(hm)*c;
      pp=p*t(p);


      wl=j(1,nq,.);
      pl=j(1,nq,.);
      prop=j(1,nq,.);
      ibin=j(1,nq,.);



      if (icontrol^=k & swtype=3 & swcov=0) then do;
        do until(icontrol>k);
          if ((swtype=4 | swtype=5) & swcov=1) then icontrol=1;
          if swtype^=0 then contr=contrast(swtype,mm,icontrol);
          gg=contr*p;
          dt=gg*t(gg);

          do i=1 to mm;
           dt[i,i]=sqrt(dt[i,i]);
          end;

          do i=1 to mm;
            do j=1 to k1;
              gg[i,j]=gg[i,j]/dt[i,i];
            end;
          end;


          reset noname;

          if (swtype=3 & swcov=0 & icontrol^=1) then do;
            run estimate(swtype,nq,jjj);

            print 'q-value  corresponding to population' icontrol ' is ' zmean;
            print mocar ' random directions used. ' 'swtype=' swtype;
            print 'Standard error of estimate is ' se;
          end;
          else do;
            if zrep<1 then jjj=-1;
            if zrep=1 then jjj=0;
            if zrep>1 then jjj=1;
            run estimate(swtype,nq,jjj);

            if (zrep<1 & swtype=zero) then do;
              print 'The estimated standard error may be very conservative';
              print 'unless the sample means are highly correlated.';
            end;
            print 'The q-value is ' zmean '. The seed is' jseed;
            print k ' populations,' conf ' confidence.';
            if ndenom>0 then print 'df for variance estimate is ' ndenom;
            else print 'variance is assumed known';
            print mocar ' random directions used. ' 'swtype=' swtype;
            print 'Standard error of estimate is ' se;
          end;
          icontrol=icontrol + 1;
        end;
      end;
      else do;
        if ((swtype=4 | swtype=5) & swcov=1) then icontrol=1;
        if swtype^=0 then contr=contrast(swtype,mm,icontrol);
        gg=contr*p;
        dt=gg*t(gg);

        do i=1 to mm;
         dt[i,i]=sqrt(dt[i,i]);
        end;

        do i=1 to mm;
          do j=1 to k1;
            gg[i,j]=gg[i,j]/dt[i,i];
          end;
        end;


        reset noname;

        if (swtype=3 & swcov=0 & icontrol^=1) then do;
          run estimate(swtype,nq,jjj);

          print 'q-value  corresponding to population' icontrol ' is ' zmean;
          print mocar ' random directions used. ' 'swtype=' swtype;
          print 'Standard error of estimate is ' se;
        end;
        else do;
          if zrep<1 then jjj=-1;
          if zrep=1 then jjj=0;
          if zrep>1 then jjj=1;
          run estimate(swtype,nq,jjj);

          if (zrep<1 & swtype=zero) then do;
            print 'The estimated standard error may be very conservative';
            print 'unless the sample means are highly correlated.';
          end;
          print 'The q-value is ' zmean '. The seed is' jseed;
          print k ' populations,' conf ' confidence.';
          if ndenom>0 then print 'df for variance estimate is ' ndenom;
          else print 'variance is assumed known';
          print mocar ' random directions used. ' 'swtype=' swtype;
          print 'Standard error of estimate is ' se;
        end;
      end;
quit;