Content-Type: text/plain; charset=US-ASCII
Content-Transfer-Encoding: 7bit
Content-MD5: xmtAukfIhMzCPOmiKcm+Jg==

Attached are the two programs for the hybrid methods.Content-Type: text/plain; charset=US-ASCII
Content-Transfer-Encoding: 7bit
Content-MD5: b9sE0rwF0tNMwUFu4cXUBQ==
Content-Disposition: inline; filename="hybrideqcov.c"
Content-Description: hybrideqcov.c
X-Content-Name: hybrideqcov.c

/*********************************************************
* C PROGRAM TO COMPARE TWO TREATMENTS ON MULTIPLE        *
* ENDPOINTS, BY COMPUTING ADJUSTED P-VALUES USING THE    *
* HYBRID METHOD OF LOGAN AND TAMHANE (2001), PUBLISHED   *
* IN THE BIOMETRICAL JOURNAL.                            *
* THIS ALGORITHM USES THE MINIMUM INDIVIDUAL P-VALUE,    *
* THE OLS P-VALUE, AND THE ALR P-VALUE, IN A CLOSED      *
* TESTING PROCEDURE.  ALL STATISTICS ARE COMPUTED        *
* UNDER THE ASSUMPTION THAT THE COVARIANCE MATRICES MAY  *
* NOT BE EQUAL.						 *
*********************************************************/


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>

/*  Command to compile in Unix:

cc -Aa hybrideqcov.c -lm

This will create an executable file called a.out.

*/

/*********************************************************
*   PARAMETERS TO BE INPUT :				 *
*     --sample sizes for each group, n_1 and n_2, where  *
*       group1 is the group expected to have a higher    *
*       response.					 *
*     --number of endpoints, cap_k			 *
*     --number of bootstrap iterations, cap_m		 *
*     --location of datafile, medata                     *
*     --location for bootstrap t and s data		 *
*							 *
*  Each observation is represented by a row in the data  *
*  file.  The first column in the data file should       *
*  indicate a 1 for the group with the higher mean, and a* 
*  2 for the other group.  The remaining cap_k columns   *
*  indicate the responses for the observation in that    *
*  row.							 *
*********************************************************/

#define n_1 54
#define n_2 57
#define cap_k 4
#define cap_m 1000
#define medata "/scr01/brl866/neqex/medata.txt"
#define boottdata "/scr01/brl866/neq/boott.dat"
#define bootsdata "/scr01/brl866/neq/boots.dat"

/* Other parameter declarations */
#define alpha 0.05
#define MAX(a,b) ((a) > (b) ? (a) :(b))
#define MIN(a,b) ((a) < (b) ? (a) :(b))
#define SQ(a) (a*a)
#define TINY 1.0e-20

/* Function declarations.  All functions follow the main
body of the program  */

void olstestctp();
void btolstestctp();
void olstest();
void olstestctpneq();
double alrtest();
void alrctp();
double tailprobt();
void normbtolsctp();
void btolstestctp();
double compute_adjp();
double tailprobf();
double bincoeff();
double betai();
double betacf();
double gammln();
void nrerror();
double ran3();

double *dvector();
int *ivector();
double **dmatrix();
int **imatrix();
double **submatrix();
void free_dvector();
void free_ivector();
void free_dmatrix();
void free_imatrix();
void free_submatrix();
double **convert_matrix();
void free_convert_matrix();
void matmult();
void scalmult();
void matadd();
void matcopy();
void ludcmp();
void lubksb();
void ludinv();
void cholesky();
void d_PrintMatrix();
void d_Transpose();
void matmean();
void matcov();
void matpoolcov();
void matsub();
void **delmatind();
void delrmatind();
void gramschmidt();

void main()
{
long *seed,temp3;
int simct,i,j,k,k2,numsig;
double **data1,**data2,*ptil,*ptilols,*ptilh,*padjalr,*padjols;
double **rawpval;
int *sigolsb,*trtmt,cnt1,cnt2;

FILE *data1file,*data2file;

/*  initialize the random number generator  */
srand(time(NULL));
temp3=-1*rand();
seed=&temp3;

printf("Hybrid Method for Equal Covariance Matrices \n\n");
data1=dmatrix(1,n_1,1,cap_k);
data2=dmatrix(1,n_2,1,cap_k);
rawpval=dmatrix(1,1,1,cap_k);
ptil=dvector(1,cap_k);
ptilh=dvector(1,cap_k);
padjalr=dvector(1,cap_k);
padjols=dvector(1,cap_k);
sigolsb=ivector(1,cap_k+1);
trtmt=ivector(1,1);

/* Read in data */

cnt1=0;
cnt2=0;
data1file=fopen(medata,"r");
  for (j=1;j<=(n_1+n_2);j++)
    {
	fscanf(data1file, "%d",&trtmt[1]);
	if (trtmt[1]==1)  {
	  cnt1+=1;
  	  for (k=1;k<=cap_k;k++)
	    {
             fscanf(data1file, "%lf",&data1[cnt1][k]);
	    }
	  }
	else {
  	  cnt2+=1;
  	  for (k=1;k<=cap_k;k++)
	    {
             fscanf(data1file, "%lf",&data2[cnt2][k]);
	    }
	  }
	
    }



/*  To print data files out

  d_PrintMatrix(data1,n_1,cap_k," %f ");
  printf("\n\n");
  d_PrintMatrix(data2,n_2,cap_k," %f ");
*/


/* Call function to apply hybrid method */

  normbtolsctp(data1,data2,n_1,n_2,cap_k,cap_m,ptil,rawpval,seed,sigolsb,ptilh);

/* Call function to apply ALR closed test procedure */

  alrctp(data1,data2,cap_k,n_1,n_2,padjalr);

/* Call function to apply OLS closed test procedure */

  olstestctp(data1,data2,cap_k,n_1,n_2,padjols);



/* Print results */

for (k=1;k<=cap_k;k++)
  {
  printf("Endpoint k=%d\n",k);
  printf("    Unadjusted p-value                          = %f\n",rawpval[1][k]);
  printf("    Adjusted p-value using max t                = %f\n",ptil[k]);
  printf("    Adjusted p-value using hybrid max t/OLS     = %f\n",ptilh[k]);
  printf("    Adjusted p-value using ALR CTP              = %f\n",padjalr[k]);
  printf("    Adjusted p-value using OLS CTP              = %f\n",padjols[k]);
  }      


free_dmatrix(data1,1,n_1,1,cap_k);
free_dmatrix(data2,1,n_2,1,cap_k);
free_dmatrix(rawpval,1,1,1,cap_k);
free_dvector(ptil,1,cap_k);
free_dvector(ptilols,1,cap_k);
free_dvector(padjalr,1,cap_k);
free_dvector(padjols,1,cap_k);
free_ivector(sigolsb,1,cap_k+1);

fclose(data1file);
fclose(data2file);

}




void normbtolsctp(double **x1,double **x2,int n1,int n2,int capk,int capm,double *origptilde,double **rawp,long *idum, int *origsigendpt,double *origptildeh)
{
int i,j,k,l,m,temp1,ranind;
double **Sstar,**S,**xbar1,**xbar2,**sorts;
double **y1,**y2,**ystarbar1,**ystarbar2,**sortxbar1,**sortxbar2;
double *t,*tstar,*sortt,*bootp,*minbootp,*ptilde,rannum,*issig,*tols;
int *rawind,*sortind,step,*sigendpt,*sigendptols,k2,df;
double **sortp,**sorty1,**sorty2,**ystar1,**ystar2,*pols,*bootpols,*minbootpols;
double **sortx1,**sortx2,*palr,*bootpalr,*sortadjph;
double tempd,w1,w2,tempd2,**mnvect1,**mnvect2,**smat1,**smat2,tempmin;
int *indendpt,sumindk;
double ran3();
void btolstestctp();
FILE *tdata,*nudata,*sdata,*s1data,*s2data,*mn1data,*mn2data;
/*
double *tempvec;
*/
rawind=ivector(1,capk);
indendpt=ivector(1,capk);
sortind=ivector(1,capk);
issig=dvector(1,capk);
Sstar=dmatrix(1,capk,1,capk);
sorts=dmatrix(1,capk,1,capk);
S=dmatrix(1,capk,1,capk);
y1=dmatrix(1,n1,1,capk);
y2=dmatrix(1,n2,1,capk);
ystar1=dmatrix(1,n1,1,capk);
ystar2=dmatrix(1,n2,1,capk);
sorty1=dmatrix(1,n1,1,capk);
sorty2=dmatrix(1,n2,1,capk);
sortx1=dmatrix(1,n1,1,capk);
sortx2=dmatrix(1,n2,1,capk);
xbar1=dmatrix(1,1,1,capk);
xbar2=dmatrix(1,1,1,capk);
ystarbar1=dmatrix(1,1,1,capk);
ystarbar2=dmatrix(1,1,1,capk);
sortxbar1=dmatrix(1,1,1,capk);
sortxbar2=dmatrix(1,1,1,capk);
t=dvector(1,capk);
tstar=dvector(1,capk);
sortt=dvector(1,capk);
bootp=dvector(1,capk);
pols=dvector(1,capk);
palr=dvector(1,capk);
bootpols=dvector(1,capk);
bootpalr=dvector(1,capk);
tols=dvector(1,capk);
minbootp=dvector(1,capk);
minbootpols=dvector(1,capk);
ptilde=dvector(1,capk);
sortadjph=dvector(1,capk);
sortp=dmatrix(1,1,1,capk);
sigendpt=ivector(1,capk);

tdata=fopen(boottdata,"w+");
sdata=fopen(bootsdata,"w+");

matmean(x1,n1,capk,xbar1);
matmean(x2,n2,capk,xbar2);
matpoolcov(x1,n1,capk,x2,n2,capk,xbar1,xbar2,S);

/* To print out sample covariance matrices or sample mean vectors 

d_PrintMatrix(S,capk,capk," %12.11f ");

d_PrintMatrix(xbar1,1,capk," %12.11f ");
d_PrintMatrix(xbar2,1,capk," %12.11f ");
*/


/*  Compute mean centered vectors to draw from */

for (j=1;j<=n1;j++)
  {
  matsub(x1+(j-1),1,capk,xbar1,1,capk,y1+(j-1));
  }
for (j=1;j<=n2;j++)
  {
  matsub(x2+(j-1),1,capk,xbar2,1,capk,y2+(j-1));
  }

/*compute rawp, order them */

for (k=1;k<=capk;k++)
  {
  t[k]=sqrt(n1*n2/(n1+n2))*(xbar1[1][k]-xbar2[1][k])/sqrt((S[k][k]));
  rawp[1][k]=tailprobt(t[k],n1+n2-2);
  rawind[k]=1;
  for (l=1;l<k;l++)
    {
    if (rawp[1][l]<rawp[1][k]) rawind[l]+=1;
    else rawind[k]+=1;
    }
  } 

  
for (k=1;k <=capk;k++)
  {
  temp1=rawind[k];
  sortind[temp1]=k;
  sortp[1][temp1]=rawp[1][k];
  sortt[temp1]=t[k];
  for (j=1;j<=n1;j++)
    {
    sorty1[j][temp1]=y1[j][k];
    sortx1[j][temp1]=x1[j][k];
    }
  for (j=1;j<=n2;j++)
    {
    sorty2[j][temp1]=y2[j][k];
    sortx2[j][temp1]=x2[j][k];
    }
  }

/* Compute ordered means, covariances, t statistics */

matmean(sortx1,n1,capk,ystarbar1);
matmean(sortx2,n2,capk,ystarbar2);
matpoolcov(sortx1,n1,capk,sortx2,n2,capk,ystarbar1,ystarbar2,sorts);

olstest(ystarbar1,ystarbar2,sorts,capk,n1,n2,tols,pols,tstar);

/* Draw capm bootstrap samples */

for (k=1;k<=capk;k++) 
  {
  issig[k]=0.0;
  }

for (m=1;m<=capm;m++)
  {
  for (j=1;j<=n1;j++)
    {
    tempd=ran3(idum);
    ranind=ceil((n1+n2)*tempd);
    if (ranind<=n1) matcopy(sorty1+ranind-1,1,capk,ystar1+j-1,1,capk);
    else matcopy(sorty2+ranind-n1-1,1,capk,ystar1+j-1,1,capk);
    }
  for (j=1;j<=n2;j++)
    {
    tempd=ran3(idum);
    ranind=ceil((n1+n2)*tempd);
    if (ranind<=n1) matcopy(sorty1+ranind-1,1,capk,ystar2+j-1,1,capk);
    else matcopy(sorty2+ranind-n1-1,1,capk,ystar2+j-1,1,capk);
    }
  
  matmean(ystar1,n1,capk,ystarbar1);
  matmean(ystar2,n2,capk,ystarbar2);
  matpoolcov(ystar1,n1,capk,ystar2,n2,capk,ystarbar1,ystarbar2,Sstar);
  olstest(ystarbar1,ystarbar2,Sstar,capk,n1,n2,tols,bootpols,tstar);

/* For each bootstrap sample, compute marginal p-value to find max t 
adjusted p-value */

  for (k=1;k<=capk;k++)
    {
    bootp[k]=tailprobt(tstar[k],n1+n2-2);


/*  Also record information for later use in a closed testing procedure */

    fprintf(tdata, "%lf ",tstar[k]);
    for (k2=1;k2<=k;k2++)
      {
      fprintf(sdata, "%lf ",Sstar[k][k2]);
      }



    if (k==1) 
      {
      minbootp[k]=bootp[k];
      }
    else 
      {
      minbootp[k]=MIN(minbootp[k-1],bootp[k]);
     }
     

    if (minbootp[k]<=sortp[1][k]) issig[k]=issig[k]+1;
    } 

  fprintf(tdata, "\n");
  fprintf(sdata,"\n");

  }

rewind(tdata);
rewind(sdata);

/*  Compute Adjusted p-values using max t approach */ 

ptilde[capk]=issig[capk]/capm;
for (k=capk-1;k>=1;k--) 
  {
  ptilde[k]=MAX(ptilde[k+1],issig[k]/capm);
  }

/*  Closed Testing Procedure for hybrid method  - calls function btolstestctp */
k=capk;
  sigendptols=ivector(1,k+1);
  for (k2=1;k2<=(k+1);k2++) sigendptols[k2]=0;

  btolstestctp(tdata,sdata,k,n1,n2,sigendptols,sortt,sorts,capk,sortadjph);

free_ivector(sigendptols,1,k+1);


/* Unsort adjusted p-values */

for (k=1;k<=capk;k++)
  {
  temp1=sortind[k];
  origptilde[temp1]=ptilde[k];
  origptildeh[temp1]=sortadjph[k];
  }

/* Free up memory*/

fclose(tdata);
fclose(sdata);


free_ivector(sigendpt,1,capk);
free_ivector(indendpt,1,capk);
free_ivector(rawind,1,capk);
free_ivector(sortind,1,capk);
free_dvector(issig,1,capk);
free_dmatrix(Sstar,1,capk,1,capk);
free_dmatrix(sorts,1,capk,1,capk);
free_dmatrix(S,1,capk,1,capk);
free_dmatrix(y1,1,n1,1,capk);
free_dmatrix(y2,1,n2,1,capk);
free_dmatrix(ystar1,1,n1,1,capk);
free_dmatrix(ystar2,1,n2,1,capk);
free_dmatrix(sorty1,1,n1,1,capk);
free_dmatrix(sorty2,1,n2,1,capk);
free_dmatrix(sortx1,1,n1,1,capk);
free_dmatrix(sortx2,1,n2,1,capk);
free_dmatrix(xbar1,1,1,1,capk);
free_dmatrix(xbar2,1,1,1,capk);
free_dmatrix(sortxbar1,1,1,1,capk);
free_dmatrix(sortxbar2,1,1,1,capk);
free_dmatrix(ystarbar1,1,1,1,capk);
free_dmatrix(ystarbar2,1,1,1,capk);
free_dvector(t,1,capk);
free_dvector(tstar,1,capk);
free_dvector(sortt,1,capk);
free_dvector(bootp,1,capk);
free_dvector(pols,1,capk);
free_dvector(palr,1,capk);
free_dvector(tols,1,capk);
free_dvector(bootpols,1,capk);
free_dvector(bootpalr,1,capk);
free_dvector(minbootp,1,capk);
free_dvector(minbootpols,1,capk);
free_dvector(ptilde,1,capk);
free_dmatrix(sortp,1,1,1,capk);

}


/* Apply Closed Test Procedure using the hybrid method,
with only the OLS test and not the ALR test */

void btolstestctp(tfile,sfile,numk,n1,n2,sigols,origt,sigmac,origk,padjh)

double *origt,**sigmac;
int numk,n1,n2, *sigols;
double *padjh;
FILE *tfile,*sfile;
{
int j2,k2,tempdf,i2,allones,allzeros,result,sumindk,*indendpt,tempi,df;
double tempnum,tempdenom,zk,tstat,tempd,tempd2;
double adjpval,*pind,w1,w2,tempnume,tempdenome;

double compute_adjp();
double tailprobt();

indendpt=ivector(1,numk);
pind=dvector(1,numk);
allones=pow(2,numk)-1;
allzeros=0;
result=pow(2,numk)-1;

/* Initial check to weed out unadjusted p-values >alpha*/

for (k2=1;k2<=numk;k2++)
  {
  pind[k2]=tailprobt(origt[k2],n1+n2-2);
  padjh[k2]=pind[k2];
  }

for (i2=(allones);i2>=1;i2--)
  {
    sumindk=0;
    for (k2=1;k2<=numk;k2++)
      {
      if (i2 & (1 << (k2-1))) indendpt[k2]=1;
      else indendpt[k2]=0;
      sumindk+=indendpt[k2];
      }
    if (sumindk>1) 
      {

      adjpval=compute_adjp(tfile,sfile,indendpt,numk,n1,n2,origt,pind,sigmac,sumindk);

      for (k2=1;k2<=numk;k2++)
        {
        if (indendpt[k2]==1) padjh[k2]=MAX(padjh[k2],adjpval);
        }
      /* to print out subset adjusted p-values 
      printf("i2=%d\n",i2);
      printf("adjpval=%f\n",adjpval);
      */  
      }
  }

free_ivector(indendpt,1,numk);
free_dvector(pind,1,numk);
}


/* Compute adjusted p-value using the hybrid test statistic
at any stage of the closed test procedure */

double compute_adjp(tfile,sfile,var_ind,kpr,n1,n2,origt,origp,origs,sumindk)
double *origt,*origp,**origs;
int kpr,n1,n2, *var_ind,sumindk;
FILE *tfile,*sfile;

{
double maxt,bootmaxt,tols,pind,pols,pmin,minp,minbootp,bootpind,bootminp;
double bootpols,*boott,**boots,sigind,padj,num,denom;
int i,j,k,m,df;
char dumstr2[640];
char *dumstr=dumstr2;

boott=dvector(1,kpr);
boots=dmatrix(1,kpr,1,kpr);
rewind(tfile);
rewind(sfile);
minp=1.0;
num=0.0;
denom=0.0;


for (i=1;i<=kpr;i++)
  {
  if (var_ind[i]==1)
    {
    
    minp=MIN(minp,origp[i]);
    num+=origt[i];
    denom+=1.0;
    for (j=1;j<i;j++)
      {
      if (var_ind[j]==1) denom=denom+2*origs[i][j]/sqrt(origs[i][i]*origs[j][j]);
      }
    }
  }
tols=num/sqrt(denom);

pols=tailprobt(tols,n1+n2-2*sumindk); /* Uses O'Brien's conservative df*/
/* pols=tailprobt(tols,100); to use Z dist */


pmin=MIN(pols,minp);

sigind=0.0;

for (m=1;m<=cap_m;m++)
  {
  for (j=1;j<=kpr;j++)
    {
    fscanf(tfile, "%lf",&boott[j]);
   for (k=1;k<=j;k++)
      {
      fscanf(sfile, "%lf",&boots[j][k]);
      }
    }
  

  bootminp=1.0;
  num=0.0;
  denom=0.0;
  for (i=1;i<=kpr;i++)
    {
    if (var_ind[i]==1)
      {
      pind=tailprobt(boott[i],n1+n2-2);
      bootminp=MIN(bootminp,pind);
      num+=boott[i];
      denom+=1.0;
      for (j=1;j<i;j++)
        {
        if (var_ind[j]==1) denom=denom+2*boots[i][j]/sqrt(boots[i][i]*boots[j][j]);
        }
      }
    }
  tols=num/sqrt(denom);
  bootpols=tailprobt(tols,n1+n2-2*sumindk); /* Uses O'Brien's conservative df */
/*  bootpols=tailprobt(tols,100); to use Z approximation */
  minbootp=MIN(bootpols,bootminp);
  if (minbootp<=pmin) sigind+=1.0;
 
  fgets(dumstr,10*8,tfile);
  fgets(dumstr,10*8*8,sfile);

  }

padj=sigind/cap_m;
free_dvector(boott,1,kpr);
free_dmatrix(boots,1,kpr,1,kpr);

return padj;
}    

/* OLS test for Equal covariance matrices */

void olstest(meanvect1,meanvect2,sigmamate,numk,n1,n2,tolsstate,pvalolse,tmarg)

double **meanvect1,**meanvect2,**sigmamate;
double *tolsstate,*pvalolse,*tmarg;
int numk,n1,n2;

{
int j2,k2,tempdf;
double tempnum,tempdenom,zk,r1,r2;
double tempnume,tempdenome;


/*
d_PrintMatrix(sigmamate,numk,numk," %12.11f ");
d_PrintMatrix(meanvect1,1,numk," %12.11f ");
d_PrintMatrix(meanvect2,1,numk," %12.11f ");
  */
tempnume=0.0;
tempdenome=0.0;
    
for (k2=1;k2<=numk;k2++)
  {
  tmarg[k2]=(meanvect1[1][k2]-meanvect2[1][k2])*sqrt(n1*n2/((n1+n2)*sigmamate[k2][k2]));
/*printf("k2=%d,tmarg=%f\n",k2,tmarg[k2]);*/
  tempnume=tempnume+tmarg[k2];
  tempdenome=tempdenome+1.0;
  for (j2=1;j2<k2;j2++)
    { 
    tempdenome=tempdenome+2*sigmamate[j2][k2]/sqrt(sigmamate[j2][j2]*sigmamate[k2][k2]);
    }
  tolsstate[k2]=tempnume/sqrt(tempdenome);
/*printf("k2=%d,tempnume=%f,tempdenome=%f,tols=%f\n",k2,tempnume,tempdenome,tolsstate[k2]);
*/
/*  pvalols[k2]=tailprobt(tolsstat[k2],n1+n2-2*numk);*/
  }
/*tols=tempnum/tempdenom;
tempdf=n1+n2-2*numk;
printf("n1=%d,n2=%d,numk=%d,tempdf=%d",n1,n2,numk,tempdf);
pvalue[numk]=tailprobt(tols,tempdf);*/
}


/* ALR test for Equal Covariance Matrices */

double alrtest(xbar1,xbar2,sigma,k,n1,n2)
double **xbar1,**xbar2,**sigma;
int k,n1,n2;
{
double **U,**A,**V,**C,**M1,**M2,**sigmainv,**tempmat1,**tempmat2,**tempmat3;
double **Y,**delta,talr,maxy,consttemp;
double pval,ftemp,fpval,temppow2,tempbin,fadj,teststat;
int i,j,nu;



U=dmatrix(1,k,1,k);
V=dmatrix(1,k,1,k);
C=dmatrix(1,k,1,k);
A=dmatrix(1,k,1,k);
M1=dmatrix(1,k,1,k);
M2=dmatrix(1,k,1,k);
sigmainv=dmatrix(1,k,1,k);
tempmat1=dmatrix(1,k,1,k);
tempmat2=dmatrix(1,k,1,k);
tempmat3=dmatrix(1,k,1,k);
Y=dmatrix(1,k,1,1);
delta=dmatrix(1,k,1,1);


/* Find A */

matcopy(sigma,k,k,tempmat1,k,k);
ludinv(tempmat1,sigmainv,k);

matcopy(sigmainv,k,k,tempmat2,k,k);
cholesky(tempmat2,C,k);

for (i=1;i<=k;i++)
  {
  for (j=1;j<=k;j++)
    {
    if (j==1) M1[i][j]=1;
    else if (i==j) M1[i][j]=1;
    else M1[i][j]=0;
    if (i==j) tempmat3[i][j]=1/sqrt(sigmainv[j][j]);
    else tempmat3[i][j]=0;
    }
  }

matcopy(C,k,k,tempmat1,k,k);
matmult(tempmat3,k,k,tempmat1,k,k,tempmat2);
ludinv(tempmat2,tempmat1,k);

for (i=1;i<=k;i++)
  {
  M2[i][1]=0.0;
  for (j=1;j<=k;j++)
    {
    M2[i][1]+=tempmat1[i][j];
    }
  for (j=2;j<=k;j++)
    {
    if (i==j) M2[i][j]=1;
    else M2[i][j]=0;
    }
  }

gramschmidt(U,M1,k);
gramschmidt(V,M2,k);

d_Transpose(V,k);
d_Transpose(C,k);
matmult(U,k,k,V,k,k,tempmat2);
matmult(tempmat2,k,k,C,k,k,A);

/* Compute test statistic */

for (i=1;i<=k;i++)
  {
  delta[i][1]=xbar1[1][i]-xbar2[1][i];
  }
matmult(A,k,k,delta,k,1,Y);
consttemp=sqrt(n1*n2/(n1+n2));
scalmult(Y,k,1,Y,consttemp);
talr=0.0;
for (i=1;i<=k;i++)
  {
  maxy=MAX(Y[i][1],0);
  talr=talr+maxy*maxy;
  }

/* COMPUTE PVAL */
teststat=talr;
pval=0.0;
nu=n1+n2-2;
for (i=1;i<=k;i++)
  {
  ftemp=teststat*(nu-k+1)/nu/i;
  fpval=tailprobf(ftemp,i,nu-k+1);

  tempbin=bincoeff(k,i);
  temppow2=pow(2,k);
  pval=pval+tempbin/temppow2*fpval;
  }
  
free_dmatrix(U,1,k,1,k);
free_dmatrix(V,1,k,1,k);
free_dmatrix(C,1,k,1,k);
free_dmatrix(A,1,k,1,k);
free_dmatrix(M1,1,k,1,k);
free_dmatrix(M2,1,k,1,k);
free_dmatrix(sigmainv,1,k,1,k);
free_dmatrix(tempmat1,1,k,1,k);
free_dmatrix(tempmat2,1,k,1,k);
free_dmatrix(tempmat3,1,k,1,k);
free_dmatrix(Y,1,k,1,1);
free_dmatrix(delta,1,k,1,1);
return pval;
}



/********************************************************
*							*
*	Incomplete Beta function and Student's T   	*
*	Distribution, due to Numerical Recipes in C	*
*							*
********************************************************/


#define NR_END 1
#define FREE_ARG char*
#define MAXIT 100
#define EPS 3.0e-7
#define FPMIN 1.0e-30

double betai(double a,double b,double x)
/*returns the incomplete beta function*/
{
  double betacf(double a,double b, double x);
  double gammln(double xx);
  void nrerror(char error_text[]);
  double bt;
  
  if (x<0.0 || x>1.0) nrerror("Bad x in routine betai");
  if (x==0.0 || x ==1.0) bt=0.0;
  else 
    bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
  if (x<((a+1.0)/(a+b+2.0)))
    return bt*betacf(a,b,x)/a;
  else 
    return 1.0-bt*betacf(b,a,1.0-x)/b;
}

double betacf(double a,double b, double x)
/*used by betai: Evaluates continued fraction for incomplete beta function
by modified Lentz's method*/
{
  void nrerror(char error_text[]);
  int m,m2;
  double aa,c,d,del,h,qab,qam,qap;
  qab=a+b;
  qap=a+1.0;
  qam=a-1.0;
  c=1.0;
  d=1.0-qab*x/qap;
  if (fabs(d) < FPMIN) d=FPMIN;
  d=1.0/d;
  h=d;
  for (m=1;m<=MAXIT;m++) {
    m2=2*m;
    aa=m*(b-m)*x/((qam+m2)*(a+m2));
    d=1.0+aa*d;
    if (fabs(d) < FPMIN) d=FPMIN;
    c=1.0+aa/c;
    if (fabs(c) < FPMIN) c=FPMIN;
    d=1.0/d;
    h *=d*c;
    aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
    d=1.0+aa*d;
    if (fabs(d)<FPMIN) d=FPMIN;
    c=1.0+aa/c;
    if (fabs(c)<FPMIN) c=FPMIN;
    d=1.0/d;
    del=d*c;
    h*=del;
    if (fabs(del-1.0)<EPS) break;
  }
  if (m>MAXIT) nrerror("a or b too big, or MAXIT too small in betacf");
  return h;
}


double gammln(double xx)
/* Returns the log of the gamma function*/
{
  double x,y,tmp,ser,tmp2;
  static double cof[6]={76.18009172947146,-86.50532032941677,
  	24.01409824083091,-1.231739572450155,0.1208650973866179e-2,
  	-0.5395239384953e-5};
  int j;
  
  y=x=xx;
  tmp=x+5.5;
  tmp-=(x+0.5)*log(tmp);
  ser=1.000000000190015;
  for (j=0;j<=5;j++) ser += cof[j]/++y;
  tmp2= -tmp+log(2.5066282746310005*ser/x);
  return tmp2;
}



void nrerror(char error_text[])
/*Numerical Recipes standard error handler */
{
  fprintf(stderr,"Numerical Recipes run-time error...\n");
  fprintf(stderr,"%s\n",error_text);
  fprintf(stderr,"...now exiting to system...\n");
  exit(1);
}

double tailprobt(double t,int dof)
{  double df,tpr,tabs,temptpr;
   df=dof+0.0;
   tabs=fabs(t);
   temptpr=betai(df/2.0,0.5,df/(df+tabs*tabs));
   if (t>0.0) tpr=temptpr/2.0;
   else tpr=1.0-temptpr/2.0;
   return tpr;
}


/* Tail probabilities of F distribution */

double tailprobf(double f,int dof1,int dof2)
{  double df1,df2,fpr;
   df1=dof1+0.0;
   df2=dof2+0.0;
   fpr=betai(df2/2.0,df1/2.0,df2/(df2+df1*f));
   return fpr;
}


/*  Binomial Coefficients, for use in ALR test */

double bincoeff(int n, int k)
{  double bnk;
  int i,num,denom;
  num=1;
  denom=1;
  for (i=n;i>(n-k);i--)
    {
    num=num*i;
    }
  for (i=1;i<=k;i++)
    {
    denom=denom*i;
    }
  bnk=num/denom;
  return bnk;
}
#undef EPS


/* Random number generator from Numerical Recipes in C 
Call with idum a negative integer to initialize */

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7 
#define RNMX (1.0-EPS)

double ran3(long *idum)

{
  int j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  double tempdoub;
  
  if (*idum <=0)  /* Initialize*/
    { 
      if (-(*idum) <1) *idum=1; /* To prevent idum=0*/
      else *idum=-(*idum);
      idum2=(*idum);
      for (j=NTAB+7;j>=0;j--)  /* Load the shuffle table,after 8 warmups*/
        {
          k=(*idum)/IQ1;
          *idum=IA1*(*idum-k*IQ1)-k*IR1;
          if (*idum<0) *idum += IM1;
          if (j < NTAB) iv[j] = *idum;
        }
      iy=iv[0];
    }
  k=(*idum)/IQ1;  /* Start here when not initializing*/
  *idum=IA1*(*idum-k*IQ1)-k*IR1;
  /*Compute idum=(IA1*idum)%IM1 without overflow by Schrage's method*/
  if (*idum < 0) *idum +=IM1;
  k=idum2/IQ2;
  idum2=IA2*(idum2-k*IQ2)-k*IR2;
  /*Compute idum2=(IA2*idum)%IM2 likewise*/
  if (idum2 < 0) idum2 +=IM2;
  j=iy/NDIV; /* will be in the range 0..NTAB-1*/
  iy=iv[j]-idum2;
  /*Here idum is shuffled, idum and idum2 are combined to generate output*/
  iv[j]= *idum;
  if (iy < 1) iy += IMM1;
  tempdoub=AM*iy;
  if (tempdoub > RNMX) tempdoub=RNMX; /*To avoid endpoint values*/
  return tempdoub;
}
#undef IM1
#undef IM2 
#undef AM 
#undef IMM1 
#undef IA1 
#undef IA2 
#undef IQ1 
#undef IQ2 
#undef IR1 
#undef IR2 
#undef NTAB 
#undef NDIV 
#undef EPS 
#undef RNMX 



/***************************************
*  Matrix Manipulation Functions Used  *
*  Throughout the Program 	       *
***************************************/

/* assigns memory for a (nh-nl+1)-dimension double precision vector */
double *dvector(nl, nh)
int nl, nh;
{
        double *v;

        v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
        if (!v) {
                printf("Run time error:\n");
                printf("allocation failure in dvector()\n");
                printf("Exit...\n");
                exit(1);
        }
        return v-nl;
}

/* assigns memory for a (nh-nl+1)-dimension integer vector */
int *ivector(nl,nh)
int nl,nh;
{
        int *v;
        v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
        if (!v) {
              printf("Run time error:\n");
              printf("allocation failure in ivector()\n");
              printf("Exit...\n");
              exit(1);
        }
        return v-nl;
}

/* assigns memory for a (nrh-nrl+1)x(nch-ncl+1) double precision matrix */
double **dmatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
        int i;
        double **m;

        m=(double **)malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
        if (!m){
                printf("Run time error:\n");
                printf("allocation failure 1 in dmatrix()\n");
                printf("Exit...\n");
                exit(1);
        }
        m-=nrl;
        for (i=nrl;i<=nrh;i++){
                m[i]=(double *)malloc((unsigned) (nch-ncl+1)*sizeof(double));
                if (!m[i]){
                        printf("Run time error:\n");
                        printf("allocation failure 2 in dmatrix()\n");
                        printf("Exit...\n");
                        exit(1);
                }
                m[i]-=ncl;
        }
        return m;
}

/* assigns memory for a (nrh-nrl+1)x(nch-ncl+1) integer matrix */
int **imatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{   
        int i;
        int **m;
   
        m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
        if (!m)
        {
                printf("Run time error:\n");
                printf("allocation failure 1 in imatrix()\n");
                printf("Exit...\n");
                exit(1);
        }
        m-=nrl;
        for (i=nrl;i<=nrh;i++)
        {
                m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
                if (!m[i])
                {
                        printf("Run time error:\n");
                        printf("allocation failure 2 in imatrix()\n");
                        printf("Exit...\n");
                        exit(1);
                } 
                m[i]-=ncl;
        }
        return m;
}
/* partition a matrix to get a submatrix */
double **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
double **a;
int oldrl,oldrh,oldcl,oldch,newrl,newcl;
{
        int i,j;
        double **m;
        m=(double **)malloc((unsigned) (oldrh-oldrl+1)*sizeof(double*));
        if (!m)
        {
                printf("Run time error:\n");
                printf("allocation failure 1 in submatrix()\n");
                printf("Exit...\n"); 
                exit(1);
        }
        m-=newrl;
        for (i=oldrl,j=newrl;i<=oldrh;i++,j++)
                m[j]=a[i]+oldcl-newcl;
        return m;         
}
/* eliminate indices of a matrix to get a submatrix */
void **delmatind(a,olddim,indexset,b,newdim)
double **a,**b;
int olddim,newdim,*indexset;
{
        int i1,i2,j1,j2;
        i2=1;
        for (i1=1;i1<=olddim;i1++)
          {
          if (indexset[i1]==1) 
            {
            b[i2][i2]=a[i1][i1];
            j2=1;
            for (j1=1;j1<i1;j1++)
              {
              if (indexset[j1]==1)
                {
                b[j2][i2]=a[j1][i1];
                b[i2][j2]=a[j1][i1];
                j2+=1;
                }
              }
            if (i2>newdim) printf("Error in dimension of submatrix\n");
            i2+=1;
            }
          }
}

/* eliminate indices of a vector to get a subvector */
void delvecind(a,olddim,indexset,b,newdim)
double *a,*b;
int olddim,newdim,*indexset;
{
        int i1,i2;
        
        i2=1;
        for (i1=1;i1<=olddim;i1++)
          {
          if (indexset[i1]==1) 
            {
            b[i2]=a[i1];
            if (i2>newdim) printf("Error in dimension of submatrix\n");
            i2+=1;
            }
          }
}

/* eliminate indices of a single row matrix to get a submatrix */
void delrmatind(a,olddim,indexset,b,newdim)
double **a,**b;
int olddim,newdim,*indexset;
{
        int i1,i2;
        
        i2=1;
        for (i1=1;i1<=olddim;i1++)
          {
          if (indexset[i1]==1) 
            {
            b[1][i2]=a[1][i1];
            if (i2>newdim) printf("Error in dimension of submatrix\n");
            i2+=1;
            }
          }
}

/* free up memory assigned to a double-precision vector */
void free_dvector(v,nl,nh)
double *v;
int nl,nh;
{
        free((char*) (v+nl));
}

/* free up memory assigned to an integer vector */
void free_ivector(v,nl,nh)
int *v,nl,nh;
{
        free((char*) (v+nl));
}

/* free up memory assigned to a double-precision vector */
void free_dmatrix(m,nrl,nrh,ncl,nch)
double **m;
int nrl,nrh,ncl,nch;
{
        int i;
        for (i=nrh; i>=nrl; i--)
                free((char*) (m[i]+ncl));
        free((char *) (m+nrl));
}

/* free up memory assigned to an integer matrix */
void free_imatrix(m,nrl,nrh,ncl,nch)
int **m;
int nrl,nrh,ncl,nch;
{   
        int i;
        for (i=nrh; i>=nrl; i--)
                free((char*) (m[i]+ncl));
        free((char *) (m+nrl));
}   

/* free up memory assigned to a submatrix */
void free_submatrix(b,nrl,nrh,ncl,nch)
double **b;
int nrl,nrh,ncl,nch;
{
        free((char*) (b+nrl));
}
/* convert a vector into a matrix */
double **convert_matrix(a,nrl,nrh,ncl,nch)
double *a;
int nrl,nrh,ncl,nch;
{
        int i,j,nrow,ncol;
        double **m;
        nrow=nrh-nrl+1;
        ncol=nch-ncl+1;
        m=(double **)malloc((unsigned) (nrow)*sizeof(double*));
        if (!m)
        {
                printf("Runtime error!\n");
                printf("allocation failure in convert_matrix()!");
                exit(1);
        }
        m-=nrl;
        for (i=1,j=nrl;i<=nrow;i++,j++)
                m[j]=a+ncol*i-ncl;
        return m;
}

/* free up memory assigned to converted matrix*/
void free_convert_matrix(b,nrl,nrh,ncl,nch)
double **b;
int nrl,nrh,ncl,nch;
{
        free((char*) (b+nrl));
}

/* perform lower-upper triangular decomposition */
void ludcmp(a,n,indx,d)
int n,*indx;
double **a,*d;
{
        int i,imax,j,k;
        double big,dum,sum,temp;
        double *vv,*dvector();
        void free_dvector();

      vv=dvector(1,n);
      *d=1.0;
      for (i=1;i<=n;i++) 
      {
                for (j=1;j<=n;j++)
              if ((temp=fabs(a[i][j])) > big) big=temp;
            if (big == 0.0) 
            {
               printf("singular matrix in routine LUDCMP()\n");
               exit(1);
            }
            vv[i]=1.0/big;
      }
      for(j=1;j<=n;j++){ 
         for (i=1;i<j;i++){   
             sum=a[i][j];
             for (k=1;k<i;k++) sum=sum-a[i][k]*a[k][j];
             a[i][j]=sum;
         }
         big=0.0;           
         for (i=j;i<=n;i++){   
             sum=a[i][j]; 
             for (k=1;k<j;k++)
                 sum=sum-a[i][k]*a[k][j];
             a[i][j]=sum;
             if  ( (dum=vv[i]*fabs(sum))  >= big) { 
                 big=dum;
                 imax=i;
             }
         }
         if ( j !=imax){
            for (k=i;k<=n;k++){ 
                dum=a[imax][k];
                a[imax][k]=a[j][k];
                a[j][k]=dum;
            }
            *d = -1.0*(*d);
            vv[imax]=vv[j];
         }
         indx[j]=imax;
         if (a[j][j] == 0.0) a[j][j]=TINY;
         if (j != n) { 
            dum=1.0/(a[j][j]);
            for (i=j+1;i<=n;i++)   a[i][j] *= dum;
         }
      }
      free_dvector(vv,1,n);
}

void lubksb(a,n,indx,b)
double **a, b[];
int n, *indx;
{
        int i,j;
        double sum;
    
        for (i=2;i<=n;i++)
        {
                sum=b[1];
                for (j=1;j<i;j++)
                        b[i]=b[i]-a[i][j]*b[j];
        }
        b[n]=b[n]/a[n][n];
        for (i=n-1; i>=1; i--)
        {
                for (j=i+1; j<=n; j++)
                        b[i]=b[i]-a[i][j]*b[j];
                b[i]=b[i]/a[i][i];
        }
}

/* using LU decomposition to invert a matrix */
void ludinv(a,y,n)
double **a, **y;
int n;
{
        double d, *col;
        int i,j,*indx;
        void free_dvector();
        double *dvector();
        int *ivector();
        void free_ivector();
    
        col=dvector(1,n);
        indx=ivector(1,n);
    
        d=1.0;
        ludcmp(a,n,indx,&d);
        for (j=1; j<=n; j++)
        {
                for (i=1; i<=n; i++)
                        col[i]=0.0;
                col[j]=1.0;
                lubksb(a,n,indx,col);
                for (i=1; i<=n; i++) 
                        y[i][j]=col[i];
        }
        free_dvector(col,1,n);
        free_ivector(indx,1,n);
}

/* performs Cholesky decomposition  */
void cholesky(a,y,n)
double **a, **y;
int n;
{
        int i,j,k;
        double sum1,sum2;
    
        for (i=1; i<=n; i++)
                for (j=n; j>i; j--)
                        y[i][j]=0.0;
    
        for (i=1; i<=n; i++)
        {
                for (j=1; j<=i; j++)
                {   
                        sum1=0.0;
                        sum2=0.0;
                        for (k=1; k<j; k++)
                        {
                                sum1=sum1+y[i][k]*y[j][k];
                                sum2=sum2+y[j][k]*y[j][k];   
                        }
                        y[i][j]=(a[i][j]-sum1)/sqrt(a[j][j]-sum2);
                }
        }
}

/* performs matrix multiplication */
void matmult(a,ra,ca,b,rb,cb,c)
double **a,**b,**c;
int ra,rb,ca,cb;
{
        int i,j,k;
    
        if (ca!=rb)
        {
                printf("Unmatched dimensions in MATMULT()!");
                exit(1);
        }
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=cb; j++)
                {
                        c[i][j]=0.0;
                        for (k=1;k<=ca;k++)
                                c[i][j]=c[i][j]+a[i][k]*b[k][j];
                }
        }
}

/* performs scalar multiplication */
void scalmult(a,ra,ca,b,c)
double **a,**b,c;
int ra,ca;
{
        int i,j,k;
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=ca; j++)
                {
                        b[i][j]=a[i][j]*c;
                }
        }
}

/* performs matrix addition */
void matadd(a,ra,ca,b,rb,cb,c)
double **a,**b,**c;
int ra,rb,ca,cb;
{
        int i,j,k;
    
        if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATadd()!");
                exit(1);
        }
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=ca; j++)
                {
                        c[i][j]=a[i][j]+b[i][j];
                }
        }
}

/* performs matrix subtraction */
void matsub(a,ra,ca,b,rb,cb,c)
double **a,**b,**c;
int ra,rb,ca,cb;
{
        int i,j,k;
    
        if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATadd()!");
                exit(1);
        }
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=ca; j++)
                {
                        c[i][j]=a[i][j]-b[i][j];
                }
        }
}

/* computes the mean vector of a matrix */
void matmean(a,ra,ca,b)
double **a,**b;
int ra,ca;
{
        int i,j,k;
     /*   if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATadd()!");
                exit(1);
        }
     */
        for (i=1;i<=ca;i++)
        {
                b[1][i]=0.0;
                for (j=1; j<=ra; j++)
                {
                        b[1][i]=a[j][i]/ra+b[1][i];
                }
        }
}

/* computes the covariance matrix of a matrix */
void matcov(a,ra,ca,b,Cov)
double **a,**b,**Cov;
int ra,ca;
{
        int i,j,k;
     /*   if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATadd()!");
                exit(1);
        }
     */
        for (i=1;i<=ca;i++)
        { 
        	for (j=1;j<=ca;j++)
        	{
        		Cov[i][j]=0.0;
        		for (k=1;k<=ra;k++)
        		{
        			Cov[i][j]+=a[k][i]*a[k][j]/(ra-1);
        		}
        		Cov[i][j]=Cov[i][j]-ra*b[1][j]*b[1][i]/(ra-1);
        	}	
        }
}

/* computes the pooled covariance matrix of two matrices */
void matpoolcov(a,ra,ca,b,rb,cb,amn,bmn,Cov)
double **a,**b,**Cov,**amn,**bmn;
int ra,ca,rb,cb;
{
        int i,j,k;
        if (ca!=cb)
        {
                printf("Unmatched dimensions in MATpoolcov()!");
                exit(1);
        }
     
        for (i=1;i<=ca;i++)
        { 
        	for (j=1;j<=ca;j++)
        	{
        		Cov[i][j]=0.0;
        		for (k=1;k<=ra;k++)
        		{
        			Cov[i][j]+=a[k][i]*a[k][j];
        		}
        		Cov[i][j]=Cov[i][j]-ra*amn[1][j]*amn[1][i];
        		for (k=1;k<=rb;k++)
        		{	
        			Cov[i][j]+=b[k][i]*b[k][j];
        		}
        		Cov[i][j]=Cov[i][j]-rb*bmn[1][j]*bmn[1][i];
        		Cov[i][j]=Cov[i][j]/(ra+rb-2);
        	}	
        }
}

/* copies a matrix into another */
void matcopy(a,ra,ca,b,rb,cb)
double **a,**b;
int ra,rb,ca,cb;
{
        int i,j,k;
    
        if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATcopy()!");
                exit(1);
        }
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=ca; j++)
                {
                        b[i][j]=a[i][j];
                }
        }
}

void d_PrintMatrix(double **M, int rows, int cols, char *fmt)
{
        int i, j;

        printf("d_PrintMatrix: %d rows and %d columns.\n", rows, cols);
        for (i=1; i<= rows; i++) {
                for (j=1; j<=cols; j++) printf(fmt, M[i][j]);
                printf("\n");
        }
}

void d_Transpose(double **A, int n) /* Transpose A, which is n x n doubles */
{
        int i, j;
        double tmp;

        for (i=1; i<=n; i++)
                for (j=(i+1); j<=n; j++) {
                        tmp = A[i][j];
                        A[i][j] = A[j][i];
                        A[j][i] = tmp;
                }
}

double vecnorm(double **mat1,double **mat2,int numcols,int col1,int col2)
{
  int i,j;
  double sum;
  sum=0.0;
  for (i=1;i<=numcols;i++)
    {
    sum=sum+mat1[i][col1]*mat2[i][col2];
    }
  return sum;
}

void gramschmidt(ortho,mat1,dim)
double **ortho,**mat1;
int dim;
{
  int i,j,k;
  double *tempvec,normnum,normdenom;

  tempvec=dvector(1,dim);
  
  for (i=1;i<=dim;i++)
    {
    ortho[i][1]=mat1[i][1];
    }
  
  for (j=2;j<=dim;j++)
    {
    for (i=1;i<=dim;i++)
      {
      tempvec[i]=mat1[i][j];
      }
    for (k=1;k<j;k++)
      {
      normnum=vecnorm(mat1,ortho,dim,j,k);
      normdenom=vecnorm(ortho,ortho,dim,k,k);
      for (i=1;i<=dim;i++)
        {
        tempvec[i]=tempvec[i]-normnum*ortho[i][k]/normdenom;
        }
      }
    for (i=1;i<=dim;i++)
      {
      ortho[i][j]=tempvec[i];
      }

    }
  for (j=1;j<=dim;j++)
    {
    normdenom=vecnorm(ortho,ortho,dim,j,j);
    for (i=1;i<=dim;i++)
      {
      ortho[i][j]=ortho[i][j]/sqrt(normdenom);
      }
    }
  free_dvector(tempvec,1,dim);
}    


void alrctp(double **datax1,double **datax2,int capk,int n1,int n2,double *padj)
/* Applies the closed testing procedure using the ALR test */
{
double teststat,**meanvect1,**meanvect2,**sigmamat1,**mnvect1,**mnvect2,**smat1;
int j2,k2,tempdf,i2,allones,allzeros,sumindk,*indendpt,tempi;
double tempnum,tempdenom,zk,tstat,tind;
double adjpval,pind;
int i;

meanvect1=dmatrix(1,1,1,capk);
meanvect2=dmatrix(1,1,1,capk);
sigmamat1=dmatrix(1,capk,1,capk);

matmean(datax1,n1,capk,meanvect1);
matmean(datax2,n2,capk,meanvect2);
matpoolcov(datax1,n1,capk,datax2,n2,capk,meanvect1,meanvect2,sigmamat1);


indendpt=ivector(1,capk);
allones=pow(2,capk)-1;
allzeros=0;

/* Compute marginal t statistics and p-values */

for (k2=1;k2<=capk;k2++)
  {
  tind=sqrt(n1*n2/(n1+n2))*(meanvect1[1][k2]-meanvect2[1][k2])/sqrt((sigmamat1[k2][k2]));
  pind=tailprobt(tind,n1+n2-2);
  padj[k2]=pind;
  }

/* Apply CTP */

for (i2=(allones);i2>=1;i2--)
  {
    sumindk=0;

/* select subset of endpoints */
    for (k2=1;k2<=capk;k2++)
      {
      if (i2 & (1 << (k2-1))) indendpt[k2]=1;
      else indendpt[k2]=0;
      sumindk+=indendpt[k2];
      }

/* Compute ALR test for subset */

    if (sumindk>1) 
      {
      mnvect1=dmatrix(1,1,1,sumindk);
      mnvect2=dmatrix(1,1,1,sumindk);
      smat1=dmatrix(1,sumindk,1,sumindk);
      delrmatind(meanvect1,capk,indendpt,mnvect1,sumindk);
      delrmatind(meanvect2,capk,indendpt,mnvect2,sumindk);
      delmatind(sigmamat1,capk,indendpt,smat1,sumindk);
      teststat=alrtest(mnvect1,mnvect2,smat1,sumindk,n1,n2);
      for (k2=1;k2<=capk;k2++)
        {
        if (indendpt[k2]==1) padj[k2]=MAX(padj[k2],teststat);
        }
      free_dmatrix(mnvect1,1,1,1,sumindk);
      free_dmatrix(mnvect2,1,1,1,sumindk);
      free_dmatrix(smat1,1,sumindk,1,sumindk);
      }
  }
  
free_ivector(indendpt,1,capk);
free_dmatrix(meanvect1,1,1,1,capk);
free_dmatrix(meanvect2,1,1,1,capk);
free_dmatrix(sigmamat1,1,capk,1,capk);
}


/* To apply OLS test in a closed test procedure */

void olstestctp(datax1,datax2,numk,n1,n2,padj)
double **datax1,**datax2,*padj;
int numk,n1,n2;
{
int j2,k2,tempdf,i2,allones,allzeros,resulte,resultpgt,sumindk,*indendpt,covind;
double tempnume,tempdenome,zk,tstat,*tmarg;
double pvalue,r1,r2;

double **meanvect1,**meanvect2,**sigmamat1,**sigmamat2,**sigmamate;
meanvect1=dmatrix(1,1,1,numk);
meanvect2=dmatrix(1,1,1,numk);
sigmamat1=dmatrix(1,numk,1,numk);
sigmamat2=dmatrix(1,numk,1,numk);
sigmamate=dmatrix(1,numk,1,numk);

indendpt=ivector(1,numk);
tmarg=dvector(1,numk);
allones=pow(2,numk)-1;
allzeros=0;
resulte=pow(2,numk)-1;
for (k2=1;k2<=numk;k2++)
  {
  padj[k2]=0.0;
  }

matmean(datax1,n1,numk,meanvect1);
matmean(datax2,n2,numk,meanvect2);
matpoolcov(datax1,n1,numk,datax2,n2,numk,meanvect1,meanvect2,sigmamate);


/* Apply CTP */

for (i2=allones;i2>=1;i2--)
  {
  sumindk=0;

/* Select subset of endpoints */

  for (k2=1;k2<=numk;k2++)
    {
    if (i2 & (1 << (k2-1))) indendpt[k2]=1;
    else indendpt[k2]=0;
    sumindk+=indendpt[k2];
    }
    
  tempnume=0.0;
  tempdenome=0.0;

/* Compute OLS test for subset of endpoints */

  for (k2=1;k2<=numk;k2++)
    {
    if (indendpt[k2]==1)
      {
      tmarg[k2]=(meanvect1[1][k2]-meanvect2[1][k2])*sqrt(n1*n2/((n1+n2)*sigmamate[k2][k2]));      
      tempnume=tempnume+tmarg[k2];
      tempdenome=tempdenome+1.0;
      for (j2=1;j2<k2;j2++)
        { 
        if (indendpt[j2]==1)  
          {
          tempdenome=tempdenome+2*sigmamate[j2][k2]/sqrt(sigmamate[j2][j2]*sigmamate[k2][k2]);
          }
        }
      }
    }
  tstat=tempnume/sqrt(tempdenome);
  tempdf=n1+n2-2*sumindk; /* Uses O'Brien's conservative df */
/*  tempdf=120;  to make it normal*/
  pvalue=tailprobt(tstat,tempdf);
  for (k2=1;k2<=numk;k2++)  
    {
    if (indendpt[k2]==1) padj[k2]=MAX(padj[k2],pvalue);
    }
  }
free_dmatrix(meanvect1,1,1,1,numk);
free_dmatrix(meanvect2,1,1,1,numk);
free_dmatrix(sigmamat1,1,numk,1,numk);
free_dmatrix(sigmamat2,1,numk,1,numk);
free_dmatrix(sigmamate,1,numk,1,numk);
free_ivector(indendpt,1,numk);
free_dvector(tmarg,1,numk);
  
  
  
}Content-Type: text/plain; charset=US-ASCII
Content-Transfer-Encoding: 7bit
Content-MD5: qWiNYDcz2NltjBWWiRZiYA==
Content-Disposition: inline; filename="hybriduneqcov.c"
Content-Description: hybriduneqcov.c
X-Content-Name: hybriduneqcov.c

/*********************************************************
* C PROGRAM TO COMPARE TWO TREATMENTS ON MULTIPLE        *
* ENDPOINTS, BY COMPUTING ADJUSTED P-VALUES USING THE    *
* HYBRID METHOD OF LOGAN AND TAMHANE (2001), PUBLISHED   *
* IN THE BIOMETRICAL JOURNAL.                            *
* THIS ALGORITHM USES THE MINIMUM INDIVIDUAL P-VALUE,    *
* THE OLS P-VALUE, AND THE ALR P-VALUE, IN A CLOSED      *
* TESTING PROCEDURE.  ALL STATISTICS ARE COMPUTED        *
* UNDER THE ASSUMPTION THAT THE COVARIANCE MATRICES MAY  *
* NOT BE EQUAL.						 *
*********************************************************/


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>

/*  Command to compile in Unix:

cc -Aa hybriduneqcov.c -lm

This will create an executable file called a.out.

*/

/*********************************************************
*   PARAMETERS TO BE INPUT :				 *
*     --sample sizes for each group, n_1 and n_2, where  *
*       group1 is the group expected to have a higher    *
*       response.					 *
*     --number of endpoints, cap_k			 *
*     --number of bootstrap iterations, cap_m		 *
*     --location of datafile, medata                     *
*     --location used for various bootstrap results      *
*							 *
*  Each observation is represented by a row in the data  *
*  file.  The first column in the data file should       *
*  indicate a 1 for the group with the higher mean, and a* 
*  2 for the other group.  The remaining cap_k columns   *
*  indicate the responses for the observation in that    *
*  row.							 *
*********************************************************/

#define n_1 54
#define n_2 57
#define cap_k 4
#define cap_m 500
#define medata "/scr01/brl866/neqex/medata.txt"
#define boottdata "/scr01/brl866/neq/boott.dat"
#define bootmn1data "/scr01/brl866/neq/bootmn1.dat"
#define bootmn2data "/scr01/brl866/neq/bootmn2.dat"
#define bootnudata "/scr01/brl866/neq/bootnu.dat"
#define bootsdata "/scr01/brl866/neq/boots.dat"
#define boots1data "/scr01/brl866/neq/boots1.dat"
#define boots2data "/scr01/brl866/neq/boots2.dat"

/* Other parameter declarations */
#define alpha 0.05
#define MAX(a,b) ((a) > (b) ? (a) :(b))
#define MIN(a,b) ((a) < (b) ? (a) :(b))
#define SQ(a) (a*a)
#define TINY 1.0e-20

/* Function declarations.  All functions follow the main
body of the program  */

void olstestctp();
void btolstestctp();
void olstestneq();
void olstestctpneq();
double alrtest1();
void alrctp();
double tailprobt();
void normbtolsctp();
void btolstestctp();
double compute_adjp();
double tailprobf();
double bincoeff();
double betai();
double betacf();
double gammln();
void nrerror();
double ran3();

double *dvector();
int *ivector();
double **dmatrix();
int **imatrix();
double **submatrix();
void free_dvector();
void free_ivector();
void free_dmatrix();
void free_imatrix();
void free_submatrix();
double **convert_matrix();
void free_convert_matrix();
void matmult();
void scalmult();
void matadd();
void matcopy();
void ludcmp();
void lubksb();
void ludinv();
void cholesky();
void d_PrintMatrix();
void d_Transpose();
void matmean();
void matcov();
void matpoolcov();
void matsub();
void **delmatind();
void delrmatind();
void gramschmidt();

void main()
{
long *seed,temp3;
int simct,i,j,k,k2,numsig;
double **data1,**data2,*ptil,*ptilols,*ptilh,*padjalr,*padjols;
double **rawpval;
int *sigolsb,*trtmt,cnt1,cnt2;

FILE *data1file,*data2file;

/*  initialize the random number generator  */
srand(time(NULL));
temp3=-1*rand();
seed=&temp3;

printf("Hybrid Method for Unequal Covariance Matrices \n\n");
data1=dmatrix(1,n_1,1,cap_k);
data2=dmatrix(1,n_2,1,cap_k);
rawpval=dmatrix(1,1,1,cap_k);
ptil=dvector(1,cap_k);
ptilh=dvector(1,cap_k);
padjalr=dvector(1,cap_k);
padjols=dvector(1,cap_k);
sigolsb=ivector(1,cap_k+1);
trtmt=ivector(1,1);

/* Read in data */

cnt1=0;
cnt2=0;
data1file=fopen(medata,"r");
  for (j=1;j<=(n_1+n_2);j++)
    {
	fscanf(data1file, "%d",&trtmt[1]);
	if (trtmt[1]==1)  {
	  cnt1+=1;
  	  for (k=1;k<=cap_k;k++)
	    {
             fscanf(data1file, "%lf",&data1[cnt1][k]);
	    }
	  }
	else {
  	  cnt2+=1;
  	  for (k=1;k<=cap_k;k++)
	    {
             fscanf(data1file, "%lf",&data2[cnt2][k]);
	    }
	  }
	
    }



/*  To print data files out

  d_PrintMatrix(data1,n_1,cap_k," %f ");
  printf("\n\n");
  d_PrintMatrix(data2,n_2,cap_k," %f ");
*/


/* Call function to apply hybrid method */

  normbtolsctp(data1,data2,n_1,n_2,cap_k,cap_m,ptil,rawpval,seed,sigolsb,ptilh);

/* Call function to apply ALR closed test procedure */

  alrctp(data1,data2,cap_k,n_1,n_2,padjalr);

/* Call function to apply OLS closed test procedure */

  olstestctpneq(data1,data2,cap_k,n_1,n_2,padjols);



/* Print results */

for (k=1;k<=cap_k;k++)
  {
  printf("Endpoint k=%d\n",k);
  printf("    Unadjusted p-value                          = %f\n",rawpval[1][k]);
  printf("    Adjusted p-value using max t                = %f\n",ptil[k]);
  printf("    Adjusted p-value using hybrid max t/OLS/ALR = %f\n",ptilh[k]);
  printf("    Adjusted p-value using ALR CTP              = %f\n",padjalr[k]);
  printf("    Adjusted p-value using OLS CTP              = %f\n",padjols[k]);
  }      


free_dmatrix(data1,1,n_1,1,cap_k);
free_dmatrix(data2,1,n_2,1,cap_k);
free_dmatrix(rawpval,1,1,1,cap_k);
free_dvector(ptil,1,cap_k);
free_dvector(ptilols,1,cap_k);
free_dvector(padjalr,1,cap_k);
free_dvector(padjols,1,cap_k);
free_ivector(sigolsb,1,cap_k+1);

fclose(data1file);
fclose(data2file);

}




void normbtolsctp(double **x1,double **x2,int n1,int n2,int capk,int capm,double *origptilde,double **rawp,long *idum, int *origsigendpt,double *origptildeh)
{
int i,j,k,l,m,temp1,ranind;
double **Sstar1,**Sstar2,**Sstare,**S1,**S2,**xbar1,**xbar2,**sorts1,**sorts2,**sortse;
double **y1,**y2,**ystarbar1,**ystarbar2,**sortxbar1,**sortxbar2;
double *t,*tstar,*sortt,*bootp,*minbootp,*ptilde,rannum,*issig,*tols;
int *rawind,*sortind,step,*sigendpt,*sigendptols,k2,df;
double **sortp,**sorty1,**sorty2,**ystar1,**ystar2,*pols,*bootpols,*minbootpols;
double **sortx1,**sortx2,*palr,*bootpalr,*sortadjph;
double tempd,w1,w2,tempd2,**mnvect1,**mnvect2,**smat1,**smat2,tempmin;
int *indendpt,sumindk;
double ran3();
void btolstestctp();
FILE *tdata,*nudata,*sdata,*s1data,*s2data,*mn1data,*mn2data;
/*
double *tempvec;
*/
rawind=ivector(1,capk);
indendpt=ivector(1,capk);
sortind=ivector(1,capk);
issig=dvector(1,capk);
Sstar1=dmatrix(1,capk,1,capk);
Sstar2=dmatrix(1,capk,1,capk);
Sstare=dmatrix(1,capk,1,capk);
sorts1=dmatrix(1,capk,1,capk);
sorts2=dmatrix(1,capk,1,capk);
sortse=dmatrix(1,capk,1,capk);
S1=dmatrix(1,capk,1,capk);
S2=dmatrix(1,capk,1,capk);
y1=dmatrix(1,n1,1,capk);
y2=dmatrix(1,n2,1,capk);
ystar1=dmatrix(1,n1,1,capk);
ystar2=dmatrix(1,n2,1,capk);
sorty1=dmatrix(1,n1,1,capk);
sorty2=dmatrix(1,n2,1,capk);
sortx1=dmatrix(1,n1,1,capk);
sortx2=dmatrix(1,n2,1,capk);
xbar1=dmatrix(1,1,1,capk);
xbar2=dmatrix(1,1,1,capk);
ystarbar1=dmatrix(1,1,1,capk);
ystarbar2=dmatrix(1,1,1,capk);
sortxbar1=dmatrix(1,1,1,capk);
sortxbar2=dmatrix(1,1,1,capk);
t=dvector(1,capk);
tstar=dvector(1,capk);
sortt=dvector(1,capk);
bootp=dvector(1,capk);
pols=dvector(1,capk);
palr=dvector(1,capk);
bootpols=dvector(1,capk);
bootpalr=dvector(1,capk);
tols=dvector(1,capk);
minbootp=dvector(1,capk);
minbootpols=dvector(1,capk);
ptilde=dvector(1,capk);
sortadjph=dvector(1,capk);
sortp=dmatrix(1,1,1,capk);
sigendpt=ivector(1,capk);

tdata=fopen(boottdata,"w+");
mn1data=fopen(bootmn1data,"w+");
mn2data=fopen(bootmn2data,"w+");
nudata=fopen(bootnudata,"w+");
sdata=fopen(bootsdata,"w+");
s1data=fopen(boots1data,"w+");
s2data=fopen(boots2data,"w+");

matmean(x1,n1,capk,xbar1);
matmean(x2,n2,capk,xbar2);
matcov(x1,n1,capk,xbar1,S1);
matcov(x2,n2,capk,xbar2,S2);

/* To print out sample covariance matrices or sample mean vectors 

d_PrintMatrix(S1,capk,capk," %12.11f ");
d_PrintMatrix(S2,capk,capk," %12.11f ");

d_PrintMatrix(xbar1,1,capk," %12.11f ");
d_PrintMatrix(xbar2,1,capk," %12.11f ");

*/

/*  Compute mean centered vectors to draw from */

for (j=1;j<=n1;j++)
  {
  matsub(x1+(j-1),1,capk,xbar1,1,capk,y1+(j-1));
  }
for (j=1;j<=n2;j++)
  {
  matsub(x2+(j-1),1,capk,xbar2,1,capk,y2+(j-1));
  }

/*compute rawp, order them */

for (k=1;k<=capk;k++)
  {
  t[k]=(xbar1[1][k]-xbar2[1][k])/sqrt((S1[k][k]/n1+S2[k][k]/n2));
  w1=S1[k][k]/n1;
  w2=S2[k][k]/n2;
  tempd2=w1+w2;
  tempd=SQ(tempd2)/(SQ(w1)/(n1-1)+SQ(w2)/(n2-1));
  df=ceil(tempd);
  rawp[1][k]=tailprobt(t[k],df);
  rawind[k]=1;
  for (l=1;l<k;l++)
    {
    if (rawp[1][l]<rawp[1][k]) rawind[l]+=1;
    else rawind[k]+=1;
    }
  } 

  
for (k=1;k <=capk;k++)
  {
  temp1=rawind[k];
  sortind[temp1]=k;
  sortp[1][temp1]=rawp[1][k];
  sortt[temp1]=t[k];
  for (j=1;j<=n1;j++)
    {
    sorty1[j][temp1]=y1[j][k];
    sortx1[j][temp1]=x1[j][k];
    }
  for (j=1;j<=n2;j++)
    {
    sorty2[j][temp1]=y2[j][k];
    sortx2[j][temp1]=x2[j][k];
    }
  }

/* Compute ordered means, covariances, t statistics */

matmean(sortx1,n1,capk,ystarbar1);
matmean(sortx2,n2,capk,ystarbar2);
matcov(sortx1,n1,capk,ystarbar1,sorts1);
matcov(sortx2,n2,capk,ystarbar2,sorts2);
matcopy(ystarbar1,1,capk,sortxbar1,1,capk);
matcopy(ystarbar2,1,capk,sortxbar2,1,capk);


olstestneq(ystarbar1,ystarbar2,sorts1,sorts2,sortse,capk,n1,n2,tols,pols,tstar);

/* Draw capm bootstrap samples */

for (k=1;k<=capk;k++) 
  {
  issig[k]=0.0;
  }

for (m=1;m<=capm;m++)
  {
  for (j=1;j<=n1;j++)
    {
    tempd=ran3(idum);
    ranind=ceil(n1*tempd);
    matcopy(sorty1+ranind-1,1,capk,ystar1+j-1,1,capk);
    }
  for (j=1;j<=n2;j++)
    {
    tempd=ran3(idum);
    ranind=ceil(n2*tempd);
    matcopy(sorty2+ranind-1,1,capk,ystar2+j-1,1,capk);
    }
  
  matmean(ystar1,n1,capk,ystarbar1);
  matmean(ystar2,n2,capk,ystarbar2);
  matcov(ystar1,n1,capk,ystarbar1,Sstar1);
  matcov(ystar2,n2,capk,ystarbar2,Sstar2);

  olstestneq(ystarbar1,ystarbar2,Sstar1,Sstar2,Sstare,capk,n1,n2,tols,bootpols,tstar);

/* For each bootstrap sample, compute marginal p-value to find max t 
adjusted p-value */

  for (k=1;k<=capk;k++)
    {
    w1=Sstar1[k][k]/n1;
    w2=Sstar2[k][k]/n2;
    tempd2=w1+w2;
    tempd=SQ(tempd2)/(SQ(w1)/(n1-1)+SQ(w2)/(n2-1));
    df=ceil(tempd);
    bootp[k]=tailprobt(tstar[k],df);


/*  Also record information for later use in a closed testing procedure */

    fprintf(tdata, "%lf ",tstar[k]);
    fprintf(mn1data, "%lf ",ystarbar1[1][k]);
    fprintf(mn2data, "%lf ",ystarbar2[1][k]);
    fprintf(nudata, "%d ",df);
    for (k2=1;k2<=k;k2++)
      {
      fprintf(sdata, "%lf ",Sstare[k][k2]);
      fprintf(s1data, "%lf ",Sstar1[k][k2]);
      fprintf(s2data, "%lf ",Sstar2[k][k2]);
      }



    if (k==1) 
      {
      minbootp[k]=bootp[k];
      }
    else 
      {
      minbootp[k]=MIN(minbootp[k-1],bootp[k]);
     }
     

    if (minbootp[k]<=sortp[1][k]) issig[k]=issig[k]+1;
    } 

  fprintf(tdata, "\n");
  fprintf(mn1data, "\n");
  fprintf(mn2data, "\n");
  fprintf(nudata, "\n");
  fprintf(sdata,"\n");
  fprintf(s1data,"\n");
  fprintf(s2data,"\n");

  }
free_dmatrix(Sstare,1,capk,1,capk);

rewind(tdata);
rewind(mn1data);
rewind(mn2data);
rewind(nudata);
rewind(sdata);
rewind(s1data);
rewind(s2data);

/*  Compute Adjusted p-values using max t approach */ 

ptilde[capk]=issig[capk]/capm;

for (k=capk-1;k>=1;k--) 
  {
  ptilde[k]=MAX(ptilde[k+1],issig[k]/capm);
  }

/*  Closed Testing Procedure for hybrid method  - calls function btolstestctp */
k=capk;
  sigendptols=ivector(1,k+1);
  for (k2=1;k2<=(k+1);k2++) sigendptols[k2]=0;

  btolstestctp(tdata,mn1data,mn2data,nudata,sdata,s1data,s2data,k,n1,n2,sigendptols,sortt,sortxbar1,sortxbar2,sorts1,sorts2,sortse,capk,sortadjph);

free_ivector(sigendptols,1,k+1);


/* Unsort adjusted p-values */

for (k=1;k<=capk;k++)
  {
  temp1=sortind[k];
  origptilde[temp1]=ptilde[k];
  origptildeh[temp1]=sortadjph[k];
  }

/* Free up memory*/

fclose(tdata);
fclose(mn1data);
fclose(mn2data);
fclose(nudata);
fclose(sdata);
fclose(s1data);
fclose(s2data);


free_ivector(sigendpt,1,capk);
free_ivector(indendpt,1,capk);
free_ivector(rawind,1,capk);
free_ivector(sortind,1,capk);
free_dvector(issig,1,capk);
free_dmatrix(Sstar1,1,capk,1,capk);
free_dmatrix(Sstar2,1,capk,1,capk);
free_dmatrix(sorts1,1,capk,1,capk);
free_dmatrix(sorts2,1,capk,1,capk);
free_dmatrix(sortse,1,capk,1,capk);
free_dmatrix(S1,1,capk,1,capk);
free_dmatrix(S2,1,capk,1,capk);
free_dmatrix(y1,1,n1,1,capk);
free_dmatrix(y2,1,n2,1,capk);
free_dmatrix(ystar1,1,n1,1,capk);
free_dmatrix(ystar2,1,n2,1,capk);
free_dmatrix(sorty1,1,n1,1,capk);
free_dmatrix(sorty2,1,n2,1,capk);
free_dmatrix(sortx1,1,n1,1,capk);
free_dmatrix(sortx2,1,n2,1,capk);
free_dmatrix(xbar1,1,1,1,capk);
free_dmatrix(xbar2,1,1,1,capk);
free_dmatrix(sortxbar1,1,1,1,capk);
free_dmatrix(sortxbar2,1,1,1,capk);
free_dmatrix(ystarbar1,1,1,1,capk);
free_dmatrix(ystarbar2,1,1,1,capk);
free_dvector(t,1,capk);
free_dvector(tstar,1,capk);
free_dvector(sortt,1,capk);
free_dvector(bootp,1,capk);
free_dvector(pols,1,capk);
free_dvector(palr,1,capk);
free_dvector(tols,1,capk);
free_dvector(bootpols,1,capk);
free_dvector(bootpalr,1,capk);
free_dvector(minbootp,1,capk);
free_dvector(minbootpols,1,capk);
free_dvector(ptilde,1,capk);
free_dmatrix(sortp,1,1,1,capk);

}


/* Apply Closed Test Procedure using the hybrid method */

void btolstestctp(tfile,mn1file,mn2file,nufile,sfile,s1file,s2file,numk,n1,n2,sigols,origt,origmn1,origmn2,origs1,origs2,sigmac,origk,padjh)

double *origt,*origmn1,*origmn2,**origs1,**origs2,**sigmac;
int numk,n1,n2, *sigols;
double *padjh;
FILE *tfile,*mn1file,*mn2file,*nufile,*sfile,*s1file,*s2file;
{
int j2,k2,tempdf,i2,allones,allzeros,result,sumindk,*indendpt,tempi,df;
double tempnum,tempdenom,zk,tstat,tempd,tempd2;
double adjpval,*pind,w1,w2,tempnume,tempdenome;

double compute_adjp();
double tailprobt();
double alrtest1();

indendpt=ivector(1,numk);
pind=dvector(1,numk);
allones=pow(2,numk)-1;
allzeros=0;
result=pow(2,numk)-1;

/* Initial check to weed out unadjusted p-values >alpha*/

for (k2=1;k2<=numk;k2++)
  {
  w1=origs1[k2][k2]/n1;
  w2=origs2[k2][k2]/n2;
  tempd2=w1+w2;
  tempd=SQ(tempd2)/(SQ(w1)/(n1-1)+SQ(w2)/(n2-1));
  df=ceil(tempd);
  pind[k2]=tailprobt(origt[k2],df);
  padjh[k2]=pind[k2];
  }

for (i2=(allones);i2>=1;i2--)
  {
    sumindk=0;
    for (k2=1;k2<=numk;k2++)
      {
      if (i2 & (1 << (k2-1))) indendpt[k2]=1;
      else indendpt[k2]=0;
      sumindk+=indendpt[k2];
      }
    if (sumindk>1) 
      {

      adjpval=compute_adjp(tfile,mn1file,mn2file,nufile,sfile,s1file,s2file,indendpt,numk,n1,n2,origt,origmn1,origmn2,pind,sigmac,origs1,origs2,sumindk);

      for (k2=1;k2<=numk;k2++)
        {
        if (indendpt[k2]==1) padjh[k2]=MAX(padjh[k2],adjpval);
        }
      /* to print out subset adjusted p-values 
      printf("i2=%d\n",i2);
      printf("adjpval=%f\n",adjpval);
      */  
      }
  }

free_ivector(indendpt,1,numk);
free_dvector(pind,1,numk);
}


/* Compute adjusted p-value using the hybrid test statistic
at any stage of the closed test procedure */

double compute_adjp(tfile,mn1file,mn2file,nufile,sfile,s1file,s2file,var_ind,kpr,n1,n2,origt,origmn1,origmn2,origp,origse,origs1,origs2,sumindk)
double *origt,*origmn1,*origmn2,*origp,**origse,**origs1,**origs2;
int kpr,n1,n2, *var_ind,sumindk;
FILE *tfile,*mn1file,*mn2file,*nufile,*sfile,*s1file,*s2file;

{
double maxt,bootmaxt,tols,pind,pols,pmin,minp,minbootp,bootpind,bootminp;
double bootpols,*boott,**bootmn1,**bootmn2,**bootse,**boots1,**boots2,sigind,padj,num,denom,w1,w2;
double **mnvect1,**mnvect2,**smat1,**smat2,palr,bootpalr;
int i,j,k,m,df,*bootnu;
char dumstr2[640];
char *dumstr=dumstr2;

boott=dvector(1,kpr);
bootmn1=dmatrix(1,1,1,kpr);
bootmn2=dmatrix(1,1,1,kpr);
bootnu=ivector(1,kpr);
bootse=dmatrix(1,kpr,1,kpr);
boots1=dmatrix(1,kpr,1,kpr);
boots2=dmatrix(1,kpr,1,kpr);
rewind(tfile);
rewind(mn1file);
rewind(mn2file);
rewind(nufile);
rewind(sfile);
rewind(s1file);
rewind(s2file);
minp=1.0;
num=0.0;
denom=0.0;


for (i=1;i<=kpr;i++)
  {
  if (var_ind[i]==1)
    {
    
    minp=MIN(minp,origp[i]);
    num+=origt[i];
    denom+=1.0;
    for (j=1;j<i;j++)
      {
      if (var_ind[j]==1) denom=denom+2*origse[i][j]/sqrt(origse[i][i]*origse[j][j]);
      }
    }
  }
tols=num/sqrt(denom);

pols=tailprobt(tols,100);

mnvect1=dmatrix(1,1,1,sumindk);
mnvect2=dmatrix(1,1,1,sumindk);
smat1=dmatrix(1,sumindk,1,sumindk);
smat2=dmatrix(1,sumindk,1,sumindk);
delrmatind(origmn1,kpr,var_ind,mnvect1,sumindk);
delrmatind(origmn2,kpr,var_ind,mnvect2,sumindk);
delmatind(origs1,kpr,var_ind,smat1,sumindk);
delmatind(origs2,kpr,var_ind,smat2,sumindk);

palr=alrtest1(mnvect1,mnvect2,smat1,smat2,sumindk,n1,n2);

free_dmatrix(mnvect1,1,1,1,sumindk);
free_dmatrix(mnvect2,1,1,1,sumindk);
free_dmatrix(smat1,1,sumindk,1,sumindk);
free_dmatrix(smat2,1,sumindk,1,sumindk);

pmin=MIN(pols,minp);
pmin=MIN(pmin,palr);

sigind=0.0;

for (m=1;m<=cap_m;m++)
  {
  for (j=1;j<=kpr;j++)
    {
    fscanf(tfile, "%lf",&boott[j]);
    fscanf(nufile, "%d",&bootnu[j]);
   for (k=1;k<=j;k++)
      {
      fscanf(sfile, "%lf",&bootse[j][k]);
      }
    }
  

  bootminp=1.0;
  num=0.0;
  denom=0.0;
  for (i=1;i<=kpr;i++)
    {
    if (var_ind[i]==1)
      {
      pind=tailprobt(boott[i],bootnu[i]);
      bootminp=MIN(bootminp,pind);
      num+=boott[i];
      denom+=1.0;
      for (j=1;j<i;j++)
        {
        if (var_ind[j]==1) denom=denom+2*bootse[i][j]/sqrt(bootse[i][i]*bootse[j][j]);
        }
      }
    }
  tols=num/sqrt(denom);
  bootpols=tailprobt(tols,100);
  minbootp=MIN(bootpols,bootminp);
  if (minbootp<=pmin) sigind+=1.0;
  else 
    {
    for (j=1;j<=kpr;j++)
      {
      fscanf(mn1file, "%lf",&bootmn1[1][j]);
      fscanf(mn2file, "%lf",&bootmn2[1][j]);

      for (k=1;k<=j;k++)
        {
        fscanf(s1file, "%lf",&boots1[j][k]);
        if (k!=j) boots1[k][j]=boots1[j][k];
        fscanf(s2file, "%lf",&boots2[j][k]);
        if (k!=j) boots2[k][j]=boots2[j][k];
        }
      }
    mnvect1=dmatrix(1,1,1,sumindk);
    mnvect2=dmatrix(1,1,1,sumindk);
    smat1=dmatrix(1,sumindk,1,sumindk);
    smat2=dmatrix(1,sumindk,1,sumindk);
    delrmatind(bootmn1,kpr,var_ind,mnvect1,sumindk);
    delrmatind(bootmn2,kpr,var_ind,mnvect2,sumindk);
    delmatind(boots1,kpr,var_ind,smat1,sumindk);
    delmatind(boots2,kpr,var_ind,smat2,sumindk);

    bootpalr=alrtest1(mnvect1,mnvect2,smat1,smat2,sumindk,n1,n2);


    /*printf("%f  ",teststat);*/
    free_dmatrix(mnvect1,1,1,1,sumindk);
    free_dmatrix(mnvect2,1,1,1,sumindk);
    free_dmatrix(smat1,1,sumindk,1,sumindk);
    free_dmatrix(smat2,1,sumindk,1,sumindk);
    
    if (bootpalr<=pmin) sigind+=1.0;   
    }
 
  fgets(dumstr,10*8,tfile);
  fgets(dumstr,10*8,mn1file);
  fgets(dumstr,10*8,mn2file);
  fgets(dumstr,10*8,nufile);
  fgets(dumstr,10*8*8,sfile);
  fgets(dumstr,10*8*8,s1file);
  fgets(dumstr,10*8*8,s2file);

  }

padj=sigind/cap_m;
free_dvector(boott,1,kpr);
free_ivector(bootnu,1,kpr);
free_dmatrix(bootse,1,kpr,1,kpr);
free_dmatrix(bootmn1,1,1,1,kpr);
free_dmatrix(bootmn2,1,1,1,kpr);
free_dmatrix(boots1,1,kpr,1,kpr);
free_dmatrix(boots2,1,kpr,1,kpr);

return padj;
}    

/* OLS test for unequal covariance matrices */

void olstestneq(meanvect1,meanvect2,sigmamat1,sigmamat2,sigmamate,numk,n1,n2,tolsstate,pvalolse,tmarg)

double **meanvect1,**meanvect2,**sigmamat1,**sigmamat2,**sigmamate;
double *tolsstate,*pvalolse,*tmarg;
int numk,n1,n2;

{
int j2,k2,tempdf;
double tempnum,tempdenom,zk,r1,r2;
double tempnume,tempdenome;


for (k2=1;k2<=numk;k2++)
  {
  for (j2=1;j2<k2;j2++)
    {
    tempnum=(sigmamat1[j2][k2]/n1+sigmamat2[j2][k2]/n2);
    tempdenom=sqrt((sigmamat1[j2][j2]/n1+sigmamat2[j2][j2]/n2)*(sigmamat1[k2][k2]/n1+sigmamat2[k2][k2]/n2));
    sigmamate[j2][k2]=tempnum/tempdenom;
    sigmamate[k2][j2]=tempnum/tempdenom;
    }
  sigmamate[k2][k2]=1.0;
  }
/*
d_PrintMatrix(sigmamate,numk,numk," %12.11f ");
d_PrintMatrix(meanvect1,1,numk," %12.11f ");
d_PrintMatrix(meanvect2,1,numk," %12.11f ");
  */
tempnume=0.0;
tempdenome=0.0;
    
for (k2=1;k2<=numk;k2++)
  {
  tmarg[k2]=(meanvect1[1][k2]-meanvect2[1][k2])/sqrt(sigmamat1[k2][k2]/n1+sigmamat2[k2][k2]/n2);
/*printf("k2=%d,tmarg=%f\n",k2,tmarg[k2]);*/
  tempnume=tempnume+tmarg[k2];
  tempdenome=tempdenome+1.0;
  for (j2=1;j2<k2;j2++)
    { 
    tempdenome=tempdenome+2*sigmamate[j2][k2]/sqrt(sigmamate[j2][j2]*sigmamate[k2][k2]);
    }
  tolsstate[k2]=tempnume/sqrt(tempdenome);
/*printf("k2=%d,tempnume=%f,tempdenome=%f,tols=%f\n",k2,tempnume,tempdenome,tolsstate[k2]);
*/
/*  pvalols[k2]=tailprobt(tolsstat[k2],n1+n2-2*numk);*/
  }
/*tols=tempnum/tempdenom;
tempdf=n1+n2-2*numk;
printf("n1=%d,n2=%d,numk=%d,tempdf=%d",n1,n2,numk,tempdf);
pvalue[numk]=tailprobt(tols,tempdf);*/
}


/* Apply ALR test */

double alrtest1(xbar1,xbar2,sigma1,sigma2,k,n1,n2)
double **xbar1,**xbar2,**sigma1,**sigma2;
int k,n1,n2;
{
double **U,**A,**V,**C,**M1,**M2,**sigma,**sigmainv,**tempmat1,**tempmat2,**tempmat3;
double **Y,**delta,talr,maxy,consttemp,**Z,tempnum,tempdenom;
double pval,ftemp,fpval,temppow2,tempbin,fadj,teststat,denom,num1,num2,nud;
double **Ssum,**Ssuminv,**tempcol,**tempcol1,**tempcol2;
int i,j,nu;


U=dmatrix(1,k,1,k);
V=dmatrix(1,k,1,k);
C=dmatrix(1,k,1,k);
A=dmatrix(1,k,1,k);
M1=dmatrix(1,k,1,k);
M2=dmatrix(1,k,1,k);
sigmainv=dmatrix(1,k,1,k);
sigma=dmatrix(1,k,1,k);
tempmat1=dmatrix(1,k,1,k);
tempmat2=dmatrix(1,k,1,k);
tempmat3=dmatrix(1,k,1,k);
Ssum=dmatrix(1,k,1,k);
Ssuminv=dmatrix(1,k,1,k);
Y=dmatrix(1,k,1,1);
Z=dmatrix(1,k,1,1);
delta=dmatrix(1,k,1,1);
tempcol=dmatrix(1,k,1,1);
tempcol1=dmatrix(1,k,1,1);
tempcol2=dmatrix(1,k,1,1);

/* Compute Z vector and Sigma_(zi,zj) */

for (i=1;i<=k;i++) 
  {
  Z[i][1]=(xbar1[1][i]-xbar2[1][i])/sqrt(sigma1[i][i]/n_1+sigma2[i][i]/n_2);
  }

for (i=1;i<=k;i++)
  {
  for (j=1;j<i;j++)
    {
    tempnum=(sigma1[j][i]/n_1+sigma2[j][i]/n_2);
    tempdenom=sqrt((sigma1[j][j]/n_1+sigma2[j][j]/n_2)*(sigma1[i][i]/n_1+sigma2[i][i]/n_2));
    sigma[j][i]=tempnum/tempdenom;
    sigma[i][j]=tempnum/tempdenom;
    }
  sigma[i][i]=1.0;
  }

/* Find A */

matcopy(sigma,k,k,tempmat1,k,k);
ludinv(tempmat1,sigmainv,k);

matcopy(sigmainv,k,k,tempmat2,k,k);
cholesky(tempmat2,C,k);

for (i=1;i<=k;i++)
  {
  for (j=1;j<=k;j++)
    {
    if (j==1) M1[i][j]=1;
    else if (i==j) M1[i][j]=1;
    else M1[i][j]=0;
    if (i==j) tempmat3[i][j]=1/sqrt(sigmainv[j][j]);
    else tempmat3[i][j]=0;
    }
  }

matcopy(C,k,k,tempmat1,k,k);
matmult(tempmat3,k,k,tempmat1,k,k,tempmat2);
ludinv(tempmat2,tempmat1,k);

for (i=1;i<=k;i++)
  {
  M2[i][1]=0.0;
  for (j=1;j<=k;j++)
    {
    M2[i][1]+=tempmat1[i][j];
    }
  for (j=2;j<=k;j++)
    {
    if (i==j) M2[i][j]=1;
    else M2[i][j]=0;
    }
  }

gramschmidt(U,M1,k);
gramschmidt(V,M2,k);

d_Transpose(V,k);
d_Transpose(C,k);
matmult(U,k,k,V,k,k,tempmat2);
matmult(tempmat2,k,k,C,k,k,A);

/* Compute test statistic */

matmult(A,k,k,Z,k,1,Y);
talr=0.0;
for (i=1;i<=k;i++)
  {
  maxy=MAX(Y[i][1],0);
  talr=talr+maxy*maxy;
  }

/*  FIND Degrees of Freedom */

for (i=1;i<=k;i++)
  {
  delta[i][1]=xbar1[1][i]-xbar2[1][i];
  }
  matadd(sigma1,k,k,sigma2,k,k,Ssum);
  ludinv(Ssum,Ssuminv,k);
  matmult(Ssuminv,k,k,delta,k,1,tempcol);
  denom=0.0;
  for (i=1;i<=k;i++)
    { 
    denom=denom+delta[i][1]*tempcol[i][1];
    }
  matmult(sigma1,k,k,tempcol,k,1,tempcol1);
  matmult(sigma2,k,k,tempcol,k,1,tempcol2);
  num1=0.0;
  num2=0.0;
  for (i=1;i<=k;i++)
    {
    num1=num1+tempcol[i][1]*tempcol1[i][1];
    num2=num2+tempcol[i][1]*tempcol2[i][1];
    }
  nud=denom*denom/(num1*num1/(n_1-1)+num2*num2/(n_2-1));
  nu=ceil(nud);
 /* nu=n_1+n_2-2;
printf("nud=%f, nu=%d",nud,nu);*/
  
/* COMPUTE PVAL */
teststat=talr;
pval=0.0;
for (i=1;i<=k;i++)
  {
  ftemp=teststat*(nu-k+1)/nu/i;
  fpval=tailprobf(ftemp,i,nu-k+1);

  tempbin=bincoeff(k,i);
  temppow2=pow(2,k);
  pval=pval+tempbin/temppow2*fpval;
  }
 
free_dmatrix(U,1,k,1,k);
free_dmatrix(V,1,k,1,k);
free_dmatrix(C,1,k,1,k);
free_dmatrix(A,1,k,1,k);
free_dmatrix(M1,1,k,1,k);
free_dmatrix(M2,1,k,1,k);
free_dmatrix(sigmainv,1,k,1,k);
free_dmatrix(sigma,1,k,1,k);
free_dmatrix(tempmat1,1,k,1,k);
free_dmatrix(tempmat2,1,k,1,k);
free_dmatrix(tempmat3,1,k,1,k);
free_dmatrix(Y,1,k,1,1);
free_dmatrix(Z,1,k,1,1);
free_dmatrix(delta,1,k,1,1);
free_dmatrix(Ssum,1,k,1,k);
free_dmatrix(Ssuminv,1,k,1,k);
free_dmatrix(tempcol,1,k,1,1);
free_dmatrix(tempcol1,1,k,1,1);
free_dmatrix(tempcol2,1,k,1,1);

return pval;
}


/********************************************************
*							*
*	Incomplete Beta function and Student's T   	*
*	Distribution, due to Numerical Recipes in C	*
*							*
********************************************************/


#define NR_END 1
#define FREE_ARG char*
#define MAXIT 100
#define EPS 3.0e-7
#define FPMIN 1.0e-30

double betai(double a,double b,double x)
/*returns the incomplete beta function*/
{
  double betacf(double a,double b, double x);
  double gammln(double xx);
  void nrerror(char error_text[]);
  double bt;
  
  if (x<0.0 || x>1.0) nrerror("Bad x in routine betai");
  if (x==0.0 || x ==1.0) bt=0.0;
  else 
    bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
  if (x<((a+1.0)/(a+b+2.0)))
    return bt*betacf(a,b,x)/a;
  else 
    return 1.0-bt*betacf(b,a,1.0-x)/b;
}

double betacf(double a,double b, double x)
/*used by betai: Evaluates continued fraction for incomplete beta function
by modified Lentz's method*/
{
  void nrerror(char error_text[]);
  int m,m2;
  double aa,c,d,del,h,qab,qam,qap;
  qab=a+b;
  qap=a+1.0;
  qam=a-1.0;
  c=1.0;
  d=1.0-qab*x/qap;
  if (fabs(d) < FPMIN) d=FPMIN;
  d=1.0/d;
  h=d;
  for (m=1;m<=MAXIT;m++) {
    m2=2*m;
    aa=m*(b-m)*x/((qam+m2)*(a+m2));
    d=1.0+aa*d;
    if (fabs(d) < FPMIN) d=FPMIN;
    c=1.0+aa/c;
    if (fabs(c) < FPMIN) c=FPMIN;
    d=1.0/d;
    h *=d*c;
    aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
    d=1.0+aa*d;
    if (fabs(d)<FPMIN) d=FPMIN;
    c=1.0+aa/c;
    if (fabs(c)<FPMIN) c=FPMIN;
    d=1.0/d;
    del=d*c;
    h*=del;
    if (fabs(del-1.0)<EPS) break;
  }
  if (m>MAXIT) nrerror("a or b too big, or MAXIT too small in betacf");
  return h;
}


double gammln(double xx)
/* Returns the log of the gamma function*/
{
  double x,y,tmp,ser,tmp2;
  static double cof[6]={76.18009172947146,-86.50532032941677,
  	24.01409824083091,-1.231739572450155,0.1208650973866179e-2,
  	-0.5395239384953e-5};
  int j;
  
  y=x=xx;
  tmp=x+5.5;
  tmp-=(x+0.5)*log(tmp);
  ser=1.000000000190015;
  for (j=0;j<=5;j++) ser += cof[j]/++y;
  tmp2= -tmp+log(2.5066282746310005*ser/x);
  return tmp2;
}



void nrerror(char error_text[])
/*Numerical Recipes standard error handler */
{
  fprintf(stderr,"Numerical Recipes run-time error...\n");
  fprintf(stderr,"%s\n",error_text);
  fprintf(stderr,"...now exiting to system...\n");
  exit(1);
}

double tailprobt(double t,int dof)
{  double df,tpr,tabs,temptpr;
   df=dof+0.0;
   tabs=fabs(t);
   temptpr=betai(df/2.0,0.5,df/(df+tabs*tabs));
   if (t>0.0) tpr=temptpr/2.0;
   else tpr=1.0-temptpr/2.0;
   return tpr;
}


/* Tail probabilities of F distribution */

double tailprobf(double f,int dof1,int dof2)
{  double df1,df2,fpr;
   df1=dof1+0.0;
   df2=dof2+0.0;
   fpr=betai(df2/2.0,df1/2.0,df2/(df2+df1*f));
   return fpr;
}


/*  Binomial Coefficients, for use in ALR test */

double bincoeff(int n, int k)
{  double bnk;
  int i,num,denom;
  num=1;
  denom=1;
  for (i=n;i>(n-k);i--)
    {
    num=num*i;
    }
  for (i=1;i<=k;i++)
    {
    denom=denom*i;
    }
  bnk=num/denom;
  return bnk;
}
#undef EPS


/* Random number generator from Numerical Recipes in C 
Call with idum a negative integer to initialize */

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7 
#define RNMX (1.0-EPS)

double ran3(long *idum)

{
  int j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  double tempdoub;
  
  if (*idum <=0)  /* Initialize*/
    { 
      if (-(*idum) <1) *idum=1; /* To prevent idum=0*/
      else *idum=-(*idum);
      idum2=(*idum);
      for (j=NTAB+7;j>=0;j--)  /* Load the shuffle table,after 8 warmups*/
        {
          k=(*idum)/IQ1;
          *idum=IA1*(*idum-k*IQ1)-k*IR1;
          if (*idum<0) *idum += IM1;
          if (j < NTAB) iv[j] = *idum;
        }
      iy=iv[0];
    }
  k=(*idum)/IQ1;  /* Start here when not initializing*/
  *idum=IA1*(*idum-k*IQ1)-k*IR1;
  /*Compute idum=(IA1*idum)%IM1 without overflow by Schrage's method*/
  if (*idum < 0) *idum +=IM1;
  k=idum2/IQ2;
  idum2=IA2*(idum2-k*IQ2)-k*IR2;
  /*Compute idum2=(IA2*idum)%IM2 likewise*/
  if (idum2 < 0) idum2 +=IM2;
  j=iy/NDIV; /* will be in the range 0..NTAB-1*/
  iy=iv[j]-idum2;
  /*Here idum is shuffled, idum and idum2 are combined to generate output*/
  iv[j]= *idum;
  if (iy < 1) iy += IMM1;
  tempdoub=AM*iy;
  if (tempdoub > RNMX) tempdoub=RNMX; /*To avoid endpoint values*/
  return tempdoub;
}
#undef IM1
#undef IM2 
#undef AM 
#undef IMM1 
#undef IA1 
#undef IA2 
#undef IQ1 
#undef IQ2 
#undef IR1 
#undef IR2 
#undef NTAB 
#undef NDIV 
#undef EPS 
#undef RNMX 



/***************************************
*  Matrix Manipulation Functions Used  *
*  Throughout the Program 	       *
***************************************/

/* assigns memory for a (nh-nl+1)-dimension double precision vector */
double *dvector(nl, nh)
int nl, nh;
{
        double *v;

        v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
        if (!v) {
                printf("Run time error:\n");
                printf("allocation failure in dvector()\n");
                printf("Exit...\n");
                exit(1);
        }
        return v-nl;
}

/* assigns memory for a (nh-nl+1)-dimension integer vector */
int *ivector(nl,nh)
int nl,nh;
{
        int *v;
        v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
        if (!v) {
              printf("Run time error:\n");
              printf("allocation failure in ivector()\n");
              printf("Exit...\n");
              exit(1);
        }
        return v-nl;
}

/* assigns memory for a (nrh-nrl+1)x(nch-ncl+1) double precision matrix */
double **dmatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
        int i;
        double **m;

        m=(double **)malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
        if (!m){
                printf("Run time error:\n");
                printf("allocation failure 1 in dmatrix()\n");
                printf("Exit...\n");
                exit(1);
        }
        m-=nrl;
        for (i=nrl;i<=nrh;i++){
                m[i]=(double *)malloc((unsigned) (nch-ncl+1)*sizeof(double));
                if (!m[i]){
                        printf("Run time error:\n");
                        printf("allocation failure 2 in dmatrix()\n");
                        printf("Exit...\n");
                        exit(1);
                }
                m[i]-=ncl;
        }
        return m;
}

/* assigns memory for a (nrh-nrl+1)x(nch-ncl+1) integer matrix */
int **imatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{   
        int i;
        int **m;
   
        m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
        if (!m)
        {
                printf("Run time error:\n");
                printf("allocation failure 1 in imatrix()\n");
                printf("Exit...\n");
                exit(1);
        }
        m-=nrl;
        for (i=nrl;i<=nrh;i++)
        {
                m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
                if (!m[i])
                {
                        printf("Run time error:\n");
                        printf("allocation failure 2 in imatrix()\n");
                        printf("Exit...\n");
                        exit(1);
                } 
                m[i]-=ncl;
        }
        return m;
}
/* partition a matrix to get a submatrix */
double **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
double **a;
int oldrl,oldrh,oldcl,oldch,newrl,newcl;
{
        int i,j;
        double **m;
        m=(double **)malloc((unsigned) (oldrh-oldrl+1)*sizeof(double*));
        if (!m)
        {
                printf("Run time error:\n");
                printf("allocation failure 1 in submatrix()\n");
                printf("Exit...\n"); 
                exit(1);
        }
        m-=newrl;
        for (i=oldrl,j=newrl;i<=oldrh;i++,j++)
                m[j]=a[i]+oldcl-newcl;
        return m;         
}
/* eliminate indices of a matrix to get a submatrix */
void **delmatind(a,olddim,indexset,b,newdim)
double **a,**b;
int olddim,newdim,*indexset;
{
        int i1,i2,j1,j2;
        i2=1;
        for (i1=1;i1<=olddim;i1++)
          {
          if (indexset[i1]==1) 
            {
            b[i2][i2]=a[i1][i1];
            j2=1;
            for (j1=1;j1<i1;j1++)
              {
              if (indexset[j1]==1)
                {
                b[j2][i2]=a[j1][i1];
                b[i2][j2]=a[j1][i1];
                j2+=1;
                }
              }
            if (i2>newdim) printf("Error in dimension of submatrix\n");
            i2+=1;
            }
          }
}

/* eliminate indices of a vector to get a subvector */
void delvecind(a,olddim,indexset,b,newdim)
double *a,*b;
int olddim,newdim,*indexset;
{
        int i1,i2;
        
        i2=1;
        for (i1=1;i1<=olddim;i1++)
          {
          if (indexset[i1]==1) 
            {
            b[i2]=a[i1];
            if (i2>newdim) printf("Error in dimension of submatrix\n");
            i2+=1;
            }
          }
}

/* eliminate indices of a single row matrix to get a submatrix */
void delrmatind(a,olddim,indexset,b,newdim)
double **a,**b;
int olddim,newdim,*indexset;
{
        int i1,i2;
        
        i2=1;
        for (i1=1;i1<=olddim;i1++)
          {
          if (indexset[i1]==1) 
            {
            b[1][i2]=a[1][i1];
            if (i2>newdim) printf("Error in dimension of submatrix\n");
            i2+=1;
            }
          }
}

/* free up memory assigned to a double-precision vector */
void free_dvector(v,nl,nh)
double *v;
int nl,nh;
{
        free((char*) (v+nl));
}

/* free up memory assigned to an integer vector */
void free_ivector(v,nl,nh)
int *v,nl,nh;
{
        free((char*) (v+nl));
}

/* free up memory assigned to a double-precision vector */
void free_dmatrix(m,nrl,nrh,ncl,nch)
double **m;
int nrl,nrh,ncl,nch;
{
        int i;
        for (i=nrh; i>=nrl; i--)
                free((char*) (m[i]+ncl));
        free((char *) (m+nrl));
}

/* free up memory assigned to an integer matrix */
void free_imatrix(m,nrl,nrh,ncl,nch)
int **m;
int nrl,nrh,ncl,nch;
{   
        int i;
        for (i=nrh; i>=nrl; i--)
                free((char*) (m[i]+ncl));
        free((char *) (m+nrl));
}   

/* free up memory assigned to a submatrix */
void free_submatrix(b,nrl,nrh,ncl,nch)
double **b;
int nrl,nrh,ncl,nch;
{
        free((char*) (b+nrl));
}
/* convert a vector into a matrix */
double **convert_matrix(a,nrl,nrh,ncl,nch)
double *a;
int nrl,nrh,ncl,nch;
{
        int i,j,nrow,ncol;
        double **m;
        nrow=nrh-nrl+1;
        ncol=nch-ncl+1;
        m=(double **)malloc((unsigned) (nrow)*sizeof(double*));
        if (!m)
        {
                printf("Runtime error!\n");
                printf("allocation failure in convert_matrix()!");
                exit(1);
        }
        m-=nrl;
        for (i=1,j=nrl;i<=nrow;i++,j++)
                m[j]=a+ncol*i-ncl;
        return m;
}

/* free up memory assigned to converted matrix*/
void free_convert_matrix(b,nrl,nrh,ncl,nch)
double **b;
int nrl,nrh,ncl,nch;
{
        free((char*) (b+nrl));
}

/* perform lower-upper triangular decomposition */
void ludcmp(a,n,indx,d)
int n,*indx;
double **a,*d;
{
        int i,imax,j,k;
        double big,dum,sum,temp;
        double *vv,*dvector();
        void free_dvector();

      vv=dvector(1,n);
      *d=1.0;
      for (i=1;i<=n;i++) 
      {
                for (j=1;j<=n;j++)
              if ((temp=fabs(a[i][j])) > big) big=temp;
            if (big == 0.0) 
            {
               printf("singular matrix in routine LUDCMP()\n");
               exit(1);
            }
            vv[i]=1.0/big;
      }
      for(j=1;j<=n;j++){ 
         for (i=1;i<j;i++){   
             sum=a[i][j];
             for (k=1;k<i;k++) sum=sum-a[i][k]*a[k][j];
             a[i][j]=sum;
         }
         big=0.0;           
         for (i=j;i<=n;i++){   
             sum=a[i][j]; 
             for (k=1;k<j;k++)
                 sum=sum-a[i][k]*a[k][j];
             a[i][j]=sum;
             if  ( (dum=vv[i]*fabs(sum))  >= big) { 
                 big=dum;
                 imax=i;
             }
         }
         if ( j !=imax){
            for (k=i;k<=n;k++){ 
                dum=a[imax][k];
                a[imax][k]=a[j][k];
                a[j][k]=dum;
            }
            *d = -1.0*(*d);
            vv[imax]=vv[j];
         }
         indx[j]=imax;
         if (a[j][j] == 0.0) a[j][j]=TINY;
         if (j != n) { 
            dum=1.0/(a[j][j]);
            for (i=j+1;i<=n;i++)   a[i][j] *= dum;
         }
      }
      free_dvector(vv,1,n);
}

void lubksb(a,n,indx,b)
double **a, b[];
int n, *indx;
{
        int i,j;
        double sum;
    
        for (i=2;i<=n;i++)
        {
                sum=b[1];
                for (j=1;j<i;j++)
                        b[i]=b[i]-a[i][j]*b[j];
        }
        b[n]=b[n]/a[n][n];
        for (i=n-1; i>=1; i--)
        {
                for (j=i+1; j<=n; j++)
                        b[i]=b[i]-a[i][j]*b[j];
                b[i]=b[i]/a[i][i];
        }
}

/* using LU decomposition to invert a matrix */
void ludinv(a,y,n)
double **a, **y;
int n;
{
        double d, *col;
        int i,j,*indx;
        void free_dvector();
        double *dvector();
        int *ivector();
        void free_ivector();
    
        col=dvector(1,n);
        indx=ivector(1,n);
    
        d=1.0;
        ludcmp(a,n,indx,&d);
        for (j=1; j<=n; j++)
        {
                for (i=1; i<=n; i++)
                        col[i]=0.0;
                col[j]=1.0;
                lubksb(a,n,indx,col);
                for (i=1; i<=n; i++) 
                        y[i][j]=col[i];
        }
        free_dvector(col,1,n);
        free_ivector(indx,1,n);
}

/* performs Cholesky decomposition  */
void cholesky(a,y,n)
double **a, **y;
int n;
{
        int i,j,k;
        double sum1,sum2;
    
        for (i=1; i<=n; i++)
                for (j=n; j>i; j--)
                        y[i][j]=0.0;
    
        for (i=1; i<=n; i++)
        {
                for (j=1; j<=i; j++)
                {   
                        sum1=0.0;
                        sum2=0.0;
                        for (k=1; k<j; k++)
                        {
                                sum1=sum1+y[i][k]*y[j][k];
                                sum2=sum2+y[j][k]*y[j][k];   
                        }
                        y[i][j]=(a[i][j]-sum1)/sqrt(a[j][j]-sum2);
                }
        }
}

/* performs matrix multiplication */
void matmult(a,ra,ca,b,rb,cb,c)
double **a,**b,**c;
int ra,rb,ca,cb;
{
        int i,j,k;
    
        if (ca!=rb)
        {
                printf("Unmatched dimensions in MATMULT()!");
                exit(1);
        }
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=cb; j++)
                {
                        c[i][j]=0.0;
                        for (k=1;k<=ca;k++)
                                c[i][j]=c[i][j]+a[i][k]*b[k][j];
                }
        }
}

/* performs scalar multiplication */
void scalmult(a,ra,ca,b,c)
double **a,**b,c;
int ra,ca;
{
        int i,j,k;
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=ca; j++)
                {
                        b[i][j]=a[i][j]*c;
                }
        }
}

/* performs matrix addition */
void matadd(a,ra,ca,b,rb,cb,c)
double **a,**b,**c;
int ra,rb,ca,cb;
{
        int i,j,k;
    
        if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATadd()!");
                exit(1);
        }
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=ca; j++)
                {
                        c[i][j]=a[i][j]+b[i][j];
                }
        }
}

/* performs matrix subtraction */
void matsub(a,ra,ca,b,rb,cb,c)
double **a,**b,**c;
int ra,rb,ca,cb;
{
        int i,j,k;
    
        if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATadd()!");
                exit(1);
        }
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=ca; j++)
                {
                        c[i][j]=a[i][j]-b[i][j];
                }
        }
}

/* computes the mean vector of a matrix */
void matmean(a,ra,ca,b)
double **a,**b;
int ra,ca;
{
        int i,j,k;
     /*   if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATadd()!");
                exit(1);
        }
     */
        for (i=1;i<=ca;i++)
        {
                b[1][i]=0.0;
                for (j=1; j<=ra; j++)
                {
                        b[1][i]=a[j][i]/ra+b[1][i];
                }
        }
}

/* computes the covariance matrix of a matrix */
void matcov(a,ra,ca,b,Cov)
double **a,**b,**Cov;
int ra,ca;
{
        int i,j,k;
     /*   if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATadd()!");
                exit(1);
        }
     */
        for (i=1;i<=ca;i++)
        { 
        	for (j=1;j<=ca;j++)
        	{
        		Cov[i][j]=0.0;
        		for (k=1;k<=ra;k++)
        		{
        			Cov[i][j]+=a[k][i]*a[k][j]/(ra-1);
        		}
        		Cov[i][j]=Cov[i][j]-ra*b[1][j]*b[1][i]/(ra-1);
        	}	
        }
}

/* computes the pooled covariance matrix of two matrices */
void matpoolcov(a,ra,ca,b,rb,cb,amn,bmn,Cov)
double **a,**b,**Cov,**amn,**bmn;
int ra,ca,rb,cb;
{
        int i,j,k;
        if (ca!=cb)
        {
                printf("Unmatched dimensions in MATpoolcov()!");
                exit(1);
        }
     
        for (i=1;i<=ca;i++)
        { 
        	for (j=1;j<=ca;j++)
        	{
        		Cov[i][j]=0.0;
        		for (k=1;k<=ra;k++)
        		{
        			Cov[i][j]+=a[k][i]*a[k][j];
        		}
        		Cov[i][j]=Cov[i][j]-ra*amn[1][j]*amn[1][i];
        		for (k=1;k<=rb;k++)
        		{	
        			Cov[i][j]+=b[k][i]*b[k][j];
        		}
        		Cov[i][j]=Cov[i][j]-rb*bmn[1][j]*bmn[1][i];
        		Cov[i][j]=Cov[i][j]/(ra+rb-2);
        	}	
        }
}

/* copies a matrix into another */
void matcopy(a,ra,ca,b,rb,cb)
double **a,**b;
int ra,rb,ca,cb;
{
        int i,j,k;
    
        if ((ra!=rb)|(ca!=cb))
        {
                printf("Unmatched dimensions in MATcopy()!");
                exit(1);
        }
    
        for (i=1;i<=ra;i++)
        {
                for (j=1; j<=ca; j++)
                {
                        b[i][j]=a[i][j];
                }
        }
}

void d_PrintMatrix(double **M, int rows, int cols, char *fmt)
{
        int i, j;

        printf("d_PrintMatrix: %d rows and %d columns.\n", rows, cols);
        for (i=1; i<= rows; i++) {
                for (j=1; j<=cols; j++) printf(fmt, M[i][j]);
                printf("\n");
        }
}

void d_Transpose(double **A, int n) /* Transpose A, which is n x n doubles */
{
        int i, j;
        double tmp;

        for (i=1; i<=n; i++)
                for (j=(i+1); j<=n; j++) {
                        tmp = A[i][j];
                        A[i][j] = A[j][i];
                        A[j][i] = tmp;
                }
}

double vecnorm(double **mat1,double **mat2,int numcols,int col1,int col2)
{
  int i,j;
  double sum;
  sum=0.0;
  for (i=1;i<=numcols;i++)
    {
    sum=sum+mat1[i][col1]*mat2[i][col2];
    }
  return sum;
}

void gramschmidt(ortho,mat1,dim)
double **ortho,**mat1;
int dim;
{
  int i,j,k;
  double *tempvec,normnum,normdenom;

  tempvec=dvector(1,dim);
  
  for (i=1;i<=dim;i++)
    {
    ortho[i][1]=mat1[i][1];
    }
  
  for (j=2;j<=dim;j++)
    {
    for (i=1;i<=dim;i++)
      {
      tempvec[i]=mat1[i][j];
      }
    for (k=1;k<j;k++)
      {
      normnum=vecnorm(mat1,ortho,dim,j,k);
      normdenom=vecnorm(ortho,ortho,dim,k,k);
      for (i=1;i<=dim;i++)
        {
        tempvec[i]=tempvec[i]-normnum*ortho[i][k]/normdenom;
        }
      }
    for (i=1;i<=dim;i++)
      {
      ortho[i][j]=tempvec[i];
      }

    }
  for (j=1;j<=dim;j++)
    {
    normdenom=vecnorm(ortho,ortho,dim,j,j);
    for (i=1;i<=dim;i++)
      {
      ortho[i][j]=ortho[i][j]/sqrt(normdenom);
      }
    }
  free_dvector(tempvec,1,dim);
}    


void alrctp(double **datax1,double **datax2,int capk,int n1,int n2,double *padj)
/* Applies the closed testing procedure using the ALR test */
{
double teststat,**meanvect1,**meanvect2,**sigmamat1,**sigmamat2,**mnvect1,**mnvect2,**smat1,**smat2;
int j2,k2,tempdf,i2,allones,allzeros,sumindk,*indendpt,tempi;
double tempnum,tempdenom,zk,tstat,tind;
double adjpval,pind;
int i;

meanvect1=dmatrix(1,1,1,capk);
meanvect2=dmatrix(1,1,1,capk);
sigmamat1=dmatrix(1,capk,1,capk);
sigmamat2=dmatrix(1,capk,1,capk);

matmean(datax1,n1,capk,meanvect1);
matmean(datax2,n2,capk,meanvect2);
matcov(datax1,n1,capk,meanvect1,sigmamat1);
matcov(datax2,n2,capk,meanvect2,sigmamat2);

indendpt=ivector(1,capk);
allones=pow(2,capk)-1;
allzeros=0;

/* Compute marginal t statistics and p-values */

for (k2=1;k2<=capk;k2++)
  {
  tind=(meanvect1[1][k2]-meanvect2[1][k2])/sqrt((sigmamat1[k2][k2]/n1+sigmamat2[k2][k2]/n2));
  pind=tailprobt(tind,n1+n2-2);
  padj[k2]=pind;
  }

/* Apply CTP */

for (i2=(allones);i2>=1;i2--)
  {
    sumindk=0;

/* select subset of endpoints */
    for (k2=1;k2<=capk;k2++)
      {
      if (i2 & (1 << (k2-1))) indendpt[k2]=1;
      else indendpt[k2]=0;
      sumindk+=indendpt[k2];
      }

/* Compute ALR test for subset */

    if (sumindk>1) 
      {
      mnvect1=dmatrix(1,1,1,sumindk);
      mnvect2=dmatrix(1,1,1,sumindk);
      smat1=dmatrix(1,sumindk,1,sumindk);
      smat2=dmatrix(1,sumindk,1,sumindk);
      delrmatind(meanvect1,capk,indendpt,mnvect1,sumindk);
      delrmatind(meanvect2,capk,indendpt,mnvect2,sumindk);
      delmatind(sigmamat1,capk,indendpt,smat1,sumindk);
      delmatind(sigmamat2,capk,indendpt,smat2,sumindk);
      teststat=alrtest1(mnvect1,mnvect2,smat1,smat2,sumindk,n1,n2);
      for (k2=1;k2<=capk;k2++)
        {
        if (indendpt[k2]==1) padj[k2]=MAX(padj[k2],teststat);
        }
      free_dmatrix(mnvect1,1,1,1,sumindk);
      free_dmatrix(mnvect2,1,1,1,sumindk);
      free_dmatrix(smat1,1,sumindk,1,sumindk);
      free_dmatrix(smat2,1,sumindk,1,sumindk);
      }
  }
  
free_ivector(indendpt,1,capk);
free_dmatrix(meanvect1,1,1,1,capk);
free_dmatrix(meanvect2,1,1,1,capk);
free_dmatrix(sigmamat1,1,capk,1,capk);
free_dmatrix(sigmamat2,1,capk,1,capk);
}


/* To apply OLS test in a closed test procedure */

void olstestctpneq(datax1,datax2,numk,n1,n2,padj)
double **datax1,**datax2,*padj;
int numk,n1,n2;
{
int j2,k2,tempdf,i2,allones,allzeros,resulte,resultpgt,sumindk,*indendpt,covind;
double tempnume,tempdenome,zk,tstat,*tmarg;
double pvalue,r1,r2;

double **meanvect1,**meanvect2,**sigmamat1,**sigmamat2,**sigmamate;
meanvect1=dmatrix(1,1,1,numk);
meanvect2=dmatrix(1,1,1,numk);
sigmamat1=dmatrix(1,numk,1,numk);
sigmamat2=dmatrix(1,numk,1,numk);
sigmamate=dmatrix(1,numk,1,numk);

indendpt=ivector(1,numk);
tmarg=dvector(1,numk);
allones=pow(2,numk)-1;
allzeros=0;
resulte=pow(2,numk)-1;
for (k2=1;k2<=numk;k2++)
  {
  padj[k2]=0.0;
  }

matmean(datax1,n1,numk,meanvect1);
matmean(datax2,n2,numk,meanvect2);
matcov(datax1,n1,numk,meanvect1,sigmamat1);
matcov(datax2,n2,numk,meanvect2,sigmamat2);

/* Compute appropriate covariance matrix */

for (k2=1;k2<=numk;k2++)
  {
  for (j2=1;j2<k2;j2++)
    {
    tempnume=(sigmamat1[j2][k2]/n1+sigmamat2[j2][k2]/n2);
    tempdenome=sqrt((sigmamat1[j2][j2]/n1+sigmamat2[j2][j2]/n2)*(sigmamat1[k2][k2]/n1+sigmamat2[k2][k2]/n2));
    sigmamate[j2][k2]=tempnume/tempdenome;
    sigmamate[k2][j2]=tempnume/tempdenome;
    }
  sigmamate[k2][k2]=1.0;
  }

/* Apply CTP */

for (i2=allones;i2>=1;i2--)
  {
  sumindk=0;

/* Select subset of endpoints */

  for (k2=1;k2<=numk;k2++)
    {
    if (i2 & (1 << (k2-1))) indendpt[k2]=1;
    else indendpt[k2]=0;
    sumindk+=indendpt[k2];
    }
    
  tempnume=0.0;
  tempdenome=0.0;

/* Compute OLS test for subset of endpoints */

  for (k2=1;k2<=numk;k2++)
    {
    if (indendpt[k2]==1)
      {
      tmarg[k2]=(meanvect1[1][k2]-meanvect2[1][k2])/sqrt(sigmamat1[k2][k2]/n1+sigmamat2[k2][k2]/n2);
      tempnume=tempnume+tmarg[k2];
      tempdenome=tempdenome+1.0;
      for (j2=1;j2<k2;j2++)
        { 
        if (indendpt[j2]==1)  
          {
          tempdenome=tempdenome+2*sigmamate[j2][k2]/sqrt(sigmamate[j2][j2]*sigmamate[k2][k2]);
          }
        }
      }
    }
  tstat=tempnume/sqrt(tempdenome);
  tempdf=n1+n2-2*sumindk; /* Uses O'Brien's conservative df */
/*  tempdf=120;  to make it normal*/
  pvalue=tailprobt(tstat,tempdf);
  for (k2=1;k2<=numk;k2++)  
    {
    if (indendpt[k2]==1) padj[k2]=MAX(padj[k2],pvalue);
    }
  }
free_dmatrix(meanvect1,1,1,1,numk);
free_dmatrix(meanvect2,1,1,1,numk);
free_dmatrix(sigmamat1,1,numk,1,numk);
free_dmatrix(sigmamat2,1,numk,1,numk);
free_dmatrix(sigmamate,1,numk,1,numk);
free_ivector(indendpt,1,numk);
free_dvector(tmarg,1,numk);
  
  
  
}