/* Order of events:

info1  swtype
info2  conf
info3  k
info4  ndenom
info5  irep

if(irep.lt.0)then
info6  givense

else
info7  mocar
  if(mocar.lt.0)then
  info8  irep
  info9  mocar


info10 swcov

if(swcov.eq.0)then
info11 vc(i,j)

if((swcov.eq.0).and.((swtype.eq.4).or.(swtype.eq.5)))then
info12 icontrol



if(swtype.eq.0)then
info13 mm
info14 contr


info15 seed

*/

%let givense=0;
%let swcov=1;
%let irep=0;
%let mocar=-1;
%let mm=0;
%let icontrol=.;




%window info1
  #10 @10 "This program calculates the constants necessary for"
  #11 @10 "obtaining simultaneous confidence intervals.For two-"
  #12 @10 "sided confidence intervals the form is"
  #13 @10 "(x(i)-x(j)) plus or minus sd(x(i)-x(j))*q/sqrt(2)"

  #15 @10 "The program has the capability of calculating q"
  #16 @10 "(or several q s for Hsu s MCB) for the following cases"
  #17 @10 "   1    Tukey"
  #18 @10 "   2    Bofinger"
  #19 @10 "   3    Hsu s MCB"
  #20 @10 "   4    Dunnett s comparisons with a control (1 sided)"
  #21 @10 "   5    Dunnett s comparisons with a control (2 sided)"
  #22 @10 "   6    Hayter"
  #23 @10 "   0    Arbitrary set of contrasts (1 and 2 sided)"
  #24 @10 "To exit from program, type -1 and press ENTER"

  #27 @10 "Enter your choice and press ENTER" +2 swtype 1 attr=underline;


%window info2
  #10 @10 "Enter desired confidence as a decimal (i.e .95) and"
  #11 @10 "press ENTER." +2 conf 5 attr=underline;


%window info3
  #10 @10 "Enter no. of populations or treatments, press enter." +2 k 2 attr=underline;

%window info4
  #10 @10 "Give df for error estimate (if the variance is assumed"
  #11 @10 "known type -1) and press enter" +2 ndenom 4 attr=underline;

%window info5
  #10 @10 "The calculation of q involves Monte Carlo methods."
  #11 @10 "You may select the number of random directions to use."
  #12 @10 "10000 is frequently sufficient, or you may specify the"
  #13 @10 "standard error of the calculated estimate for q"

  #15 @10 "To specify the standard error type -1 and press enter"
  #16 @10 "For no. of random directions, type 0, press enter."
  #18 @10 "Enter your choice and press ENTER" +2 irep 3 attr=underline;


%window info6
  #10 @10 "The smaller the standard error, the longer the running"
  #11 @10 "time of the program will be. Some choices are .01,"
  #12 @10 ".005, or .001. Type choice, press ENTER" +2 givense 7 attr=underline;


%window info7
  #10 @10 "Type number of random directions you wish to use and"
  #11 @10 "press enter. (The standard error of your estimate will"
  #12 @10 "be calculated using an empirical formula.)"

  #14 @10 "An alternate procedure is to make n independent esti-"
  #15 @10 "mates of q, each with nn random directions. The stand-"
  #16 @10 "ard error of estimated q is calculated using the n"
  #17 @10 "estimates. The alternate procedure is especially use-"
  #18 @10 "ful when the estimates of treatment means are"
  #19 @10 "correlated and an accurate formula for the standard"
  #20 @10 "error of the estimate of q is not available. "

  #22 @10 "If the alternate procedure is desired, type -1 and"
  #23 @10 "press enter" +2 mocar 7 attr=underline;


%window info8
  #10 @10 "Give n and press enter." +2 irep 7 attr=underline;

%window info9
  #10 @10 "Give nn and press enter." +2 mocar 7 attr=underline;



%window info10
  #10 @10 "The program will calculate the constant q for an"
  #11 @10 "arbitrary variance-covariance matrix (correlated case)"
  #12 @10 "or the uncorrelated case."

  #14 @10 "        0        correlated or unequal var case"
  #15 @10 "        1        uncorrelated equal variance case"
  #16 @10 "Enter 0 or 1 and press enter" +2 swcov 2 attr=underline;



%window info11
  #10 @10 "Input the lower triangular portion of the variance covariance matrix."
  #12 @10 "Press enter after each value" +2 vc_i 7 attr=underline;


%macro mvar;
  %display info10;
  %if (&swcov=0) %then %do;
    %do i=1 %to &k;
      %do j=1 %to &i;
         %display info11;
         vc[&i,&j]=&vc_i;
      %end;
    %end;
  %end;
  %else %put 'hallo';
%mend;



%window info12
  #10 @10 "Type number of the control population and press enter." +2 icontrol 2 attr=underline;


%window info13
  #10 @10 "A confidence interval of the form (g,+infinity) for"
  #11 @10 "-x2+x4 is input as 0 -1 0 +1  (one line)"
  #12 @10 "A confidence interval of the form (-infinity,g) for"
  #14 @10 "-x2+x4 is input as 0 +1 0 -1  (one line)"
  #14 @10 "A confidence interval of the form (-g,g) for -x2+x4"
  #15 @10 "requires both  (two lines)"
  #16 @10 "Type no. of lines required for contrasts,"
  #27 @10 "press enter" +2 mm 2 attr=underline;


%window info14
  #10 @10 "Type the contrasts, pressing enter after each value" +2 contr_i 7 attr=underline;

%macro mcontr;
  %if (&swtype=0) %then %do;
    %do i=1 %to &mm;
      %do j=1 %to &k;
         %display info14;
         contr[&i,&j]=&contr_i;
      %end;
    %end;
  %end;
  %else %put 'hallo';
%mend;


%window info15
  #10 @10 "Type seed for random no generator, press enter" +2 seed 10 attr=underline;


%macro main;
  %display info1;
  %display info2;
  %display info3;
  %display info4;
  %display info5;
  %if (&irep<0) %then %display info6;
  %else %do;
    %display info7;
    %if (&mocar<0) %then %do;
      %display info8;
      %display info9;
    %end;
  %end;
  %if (&swtype=0) %then %display info13;
  %display info15;
%mend;


%macro icontr;
  %if (&swcov=0) %then %do;
    %if ((&swtype=4) | (&swtype=5)) %then %display info12;
  %end;
%mend;



%main;


proc iml symsize=250000;

  swtype   = &swtype;
  conf     = &conf;
  k        = &k;
  ndenom   = &ndenom;
  irep     = &irep;
  mocar    = &mocar;
  givense  = &givense;


  vc=j(k,k,.);
  %mvar;
  %icontr;
  swcov   = &swcov;
  icontrol = &icontrol;
  do i=1 to k-1;
    do j=i+1 to k;
      vc[i,j]=vc[j,i];
    end;
  end;


  if (swtype=0) then do;
    mm = &mm;
    contr=j(mm,k,.);
  end;
  %mcontr;

  seed = &seed;




/*
  swcov  = 0     ;
  swtype = 3     ;
  conf   = .9    ;
  givense= 999   ;
  k      = 4     ;
  seed   = 8763  ;
  ndenom = 25    ;
  mocar  = 1000  ;
  irep   = 100   ;
  mm     = 99    ;
*/

  nq     = 25    ;


  /* Optional Parameters

  vc={76.807 36.148 -20.366 7.2384,
      36.148 27.373 15.176 21.198,
      -20.366 15.176 64.583 40.605,
      7.2384 21.198 40.605 31.239};
  contr=.;
  icontrol=.;    */







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

      /* 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(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;
print contr;
        return(contr);
      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) 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^=0)             then contr=j(mm,k,.);


      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;