/********

Version: July 26th.

Adding the parameter SINGLE (0 remove single samples, 1 leave untouched)
*********/


/************************************************************************************
 rpm.c      
 a program to read in data to be used in the RPM method
*************************************************************************************/

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <time.h>

#define PATHSIZE 100
#define MAXLABELS 100
#define LINESIZE 100000

#define ADR(t,x) ( t=x, &t )

// used to permute the traits
struct array {
  double trait;
  int perm;
};

struct space{
  char str[20]; //store a string
};

// used for the first output table 
struct table1{
  int group;
  char *label;
  int num;
};

typedef struct space   space;
typedef struct table1  table1;


void read_data(char pfile[], char ifile[], int *num_of_rows, int *tNum, int *traitType, int *cNum, 
               int *gNum, int *MissValue, int *permNum, long int *seed, int *zerorsq, float *pseq100,
	       float *pseq1000, double *alpha, int *ncomb, int *single, int *covGroups, 
               int **TraitTest, int *tCounter, int **CovTest[], int **cCounter, int **ID, 
               double **TD[], space **CG[]);


void make_label(double **trait[], int **N, char **label[], int *nlabels, int **indirect, int **invalid, 
		int nrows, double **TD, int Tcol, space **CG, int CovCols[], int nCov, char *sMissValue );


double run_rpm(double **trait, int *N, int nlabels, 
	       int *indirect, int *invalid, double alpha, int single, double **sum, double **sumsq);

int choose_merge(int nlabels, int tnlabels, double  *sum, double *sumsq, int N[], 
                 int invalid[], int *loc1, int *loc2, double alpha, int single);

void get_output(table1 **outtable, int indirect[], char *label[], int nlabels, int N[]);

double rsq(int nlabels, double  *sum, double *sumsq, int N[], int invalid[], int single);

void perm_trait(int num_of_rows, double **TD[], int Tcol, double **trait[], int *N, int nlabels);



/**************************************************************************************
 method: main                                                                       
 desc:   The main program. This method reads the text files and run the rpm program
 param:  none 
 return: 0 if exits without errors, else 1
 author: Tsvika Klein                                                               
 date:   03/25/04                                                              
 **************************************************************************************/

int main(int argc, char *argv[])
{

  int *indirect; // indicates the cells merged together. 
  int *invalid;  // indicates cells that are not active anymore (after merge)
  double **trait; //holds all the traits val for each unique label
  int *N; //number of traits(samples) for each unique label
  int Norig[MAXLABELS]; //number of traits(samples) for each unique label
  char **label; // holds all the labels

  double *sum;   // holds the sum of the traits for each unique label
  double *sumsq; // holds the sum of the traits squared for each unique label
  double var; // holds the variance of a specific group before printing to .table

  int nlabels=0; //number of unique labels

  int **outtable2; // holds the group assignment for table2

  char *s;//temp string for a new label
  char sMissValue[30];//temp string for a new label
  
  table1 *outtable; // holds the info for table1
  int i,p,t,l1,l2,l3,g,n,x,y,ind; 
  char *param = "param"; // default param name
  char *input = "data";  // default input file name
  char *output = "data.out"; // default output file name
  char output_table[200]; // store the outtable file name
  char output_table2[200]; // store the outtable file name
  FILE *ofp, *ofs; // output file point and seed file pointer
  FILE **ofp2; // another output file pointer
  FILE *ofp_table, *ofp_table2; // file pointers for table1 & table2


  int *CovCols; // the two covariates to be tested

  int num_of_rows = 0; // store total number of lines in the data
  int tNum,  // the number of traits in a line of data 
      traitType = 1,  // the quantative and qualitive representation
      cNum,  // the number of covariates in a line of data
      gNum,  // the number of genotypes in a line of data
      MissValue = -777,  // the number used to represent the missing value in the input data
      permNum = 1000,  // the number of permutations
      zerorsq = 1, //  0 - run perm for all, 1 - skip permutations if original rsq is 0   
      numOrigRsq = 0, // number of rsq to be stored in an array
      ncomb = 2, // number of combinations of loci to look at
      index = 0,
      numMissing=0,  // number of missing individuals from the analysis 
      single=0;

  float pseq100=0.1; // sequential testing use Pval<0.1 after 100 permutations
  float pseq1000=0.01; // sequential testing use Pval<0.01 after 1000 permutations
  double alpha = 0.99; // alpha to be used in tuckey

  long int seed = 0; // a seed value.


  int *ID;           // ID array
  double **TD;  // double type array for "Traits" if it is quantitive 
  double **tTD;
  struct space **CG; //store "Covariates" and "Geno Loci" 

  int *TraitTest;  // traits to use for the test
  int **CovTest;   // covariates to use for the test
  int covGroups;   // determines the number of covariate groups

  double origrsq; // the rsq values for all the combinatios
  double pval; // the p-values for all the combinations
  double permrsq; // the rsq of the permutation
 
  int tCounter, *cCounter; // number of traits and covariates.

  // read arguments in case the user specified them
  for (i=0;i<argc;i++){

    if (strcmp(argv[i],"-input") == 0){
      input = (char*)malloc( (strlen(argv[i+1])+1) * sizeof(char) );
      strcpy(input, argv[i+1]);
    }

    if (strcmp(argv[i],"-output") == 0){
      output = (char*)malloc( (strlen(argv[i+1])+1) * sizeof(char) );
      strcpy(output, argv[i+1]);
    }

    if (strcmp(argv[i],"-param") == 0){
      param = (char*)malloc( (strlen(argv[i+1])+1) * sizeof(char) );
      strcpy(param, argv[i+1]);
    }

    if (strcmp(argv[i],"-help") == 0){
      printf("You can specify parameters to the rpm as follows:\n");
      printf("rpmS.o -input [filename] -output [filename] -param [filename]\n");
      return 0;
    }
  }
  
  // hold the names for the tables
  strcpy(output_table, output);
  strcat(output_table, ".table");

  strcpy(output_table2, output);
  strcat(output_table2, ".table2");

  ofp_table = fopen(output_table,"a");
  ofp_table2 = fopen(output_table2,"a");

  read_data(param, input, &num_of_rows, &tNum, &traitType, &cNum, &gNum, &MissValue, &permNum, &seed, 
            &zerorsq, &pseq100, &pseq1000, &alpha, &ncomb, &single, &covGroups, &TraitTest, &tCounter, 
	    &CovTest, &cCounter, &ID, &TD, &CG); 

  sprintf(sMissValue,"%d",MissValue);

  CovCols = (int*) malloc(ncomb * sizeof(int));

  ofp = fopen(output,"a");
  
  // if the seed is 0 generate a new one
  if (seed == 0)
    seed = time(NULL);
  srandom(seed);

  // allocate memory for a label 
  s = (char*) malloc(((cNum+gNum)*10+(cNum+gNum)-1) * sizeof(char));

  // allocate memory for table2 and for a temporary holder of the traits
  outtable2 = (int **) calloc(num_of_rows, sizeof(int*));
  tTD = (double **) malloc(num_of_rows * sizeof(double*)); 
  for (n=0;n<num_of_rows;n++){
    (tTD)[n] = (double*) malloc(tNum * sizeof(double));
    if (ncomb == 3)
      (outtable2)[n] = (int*) calloc( (tNum * (cNum+gNum) * (cNum+gNum-1)*(cNum+gNum-2)  / 6),  sizeof(int)); 
    if (ncomb == 2)
      (outtable2)[n] = (int*) calloc( (tNum * (cNum+gNum) * (cNum+gNum-1) / 2),  sizeof(int));
    for (t=0;t<tNum;t++){ 
      tTD[n][t] = TD[n][t];
    }
  }
 
  if (ncomb == 3){ // if we look at 3 loci at a time
    ind=0;
    for (t=0;t<tCounter;t++){
      for (g=0;g<covGroups;g++){
	for (l1=0;l1<(cCounter[g]-2);l1++){
	  for (l2=(l1+1);l2<cCounter[g]-1;l2++){
	    for (l3=(l2+1);l3<cCounter[g];l3++){
	      CovCols[0] = CovTest[g][l1]; // the two cov
	      CovCols[1] = CovTest[g][l2]; // to be tested
	      CovCols[2] = CovTest[g][l3];

	      make_label(&trait, &N, &label, &nlabels, &indirect, &invalid, num_of_rows, TD, TraitTest[t], 
			 CG, CovCols, 3,sMissValue);


	      for (i=0;i<nlabels;i++)
		Norig[i] = N[i]; // an array with the original information 

	      origrsq = run_rpm(trait, N, nlabels, indirect, invalid, alpha, single, &sum, &sumsq);

	      
	      get_output(&outtable, indirect, label, nlabels, Norig);
	      
	  
	      for (i=0;i<num_of_rows;i++){
		for (n=0;n<nlabels;n++){
		  sprintf(s,"%s", CG[i][CovCols[0]].str); 
		  sprintf(s, "%s_%s", s, CG[i][CovCols[1]].str);
		  sprintf(s, "%s_%s", s, CG[i][CovCols[2]].str);
		  
		  if (strcmp(outtable[n].label,s) == 0){
		    outtable2[i][ind] = outtable[n].group;
		  }
		}
	      }
	      ind++;
	      
	      // print to file the output and table1
	      fprintf(ofp, "%d\t%d\t%d\t%d\t0\t%ld\t%f\n", 
		      TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, CovCols[2]+1, seed, origrsq);
	      
	      fprintf(ofp_table, "T=%d\tL1=%d\tL2=%d\tL3=%d\n", 
		      TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, CovCols[2]+1);
	      numMissing = num_of_rows;
	      for (i=0, n=0;i<nlabels;i++) {
		if (invalid[i] == 0){
		  n++;
		  numMissing-=N[i];
		  var = (sumsq[i]-(pow(sum[i],2)/N[i]))/(N[i]-1);
		  if (single==0 && N[i]==1)
		    fprintf(ofp_table, "G%d: mean=%f  sd=%f  n=%d  (ignored in rsq)\n",
			    n,sum[i]/N[i],sqrt(var),N[i]);
		  else 
		    fprintf(ofp_table, "G%d: mean=%f  sd=%f  n=%d  (included in rsq)\n",
			    n,sum[i]/N[i],sqrt(var),N[i]);
		}
	      }
	      fprintf(ofp_table, "Missing=%d\n",numMissing);
	
	      pval = 0;
	      if (origrsq != 0 || zerorsq == 0){

		// first 100 permutations or less
		for (p=1;p<=permNum && p<=100;p++){
		
		  for (i=0;i<nlabels;i++){
		    N[i] = Norig[i];
		    indirect[i] = i;
		    invalid[i] = 0;
		  }
		  // premute the trait and calculate the permuted rsq after running rpm
		  perm_trait(num_of_rows, &tTD, TraitTest[t], &trait, N, nlabels); 
	    
		  permrsq = run_rpm(trait, N, nlabels, indirect, invalid, alpha, single, &sum, &sumsq);
	   
		  if (origrsq <= permrsq)
		    pval++; 

		  fprintf(ofp, "%d\t%d\t%d\t%d\t%d\t%ld\t%f\n", 
			  TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, CovCols[2]+1, p, seed, permrsq);
		}

		// check pval after 100 permutations, if satisfies the condition
		// run up to 1000 permutations
		if ( (pval/100)<=pseq100 && permNum>=100 ){

		  for (p=101;p<=permNum && p<=1000;p++){

		    for (i=0;i<nlabels;i++){
		      N[i] = Norig[i];
		      indirect[i] = i;
		      invalid[i] = 0;
		    }
		    // premute the trait and calculate the permuted rsq after running rpm
		    perm_trait(num_of_rows, &tTD, TraitTest[t], &trait, N, nlabels); 
	    
		    permrsq = run_rpm(trait, N, nlabels, indirect, invalid, alpha, single, &sum, &sumsq);
		    
		    if (origrsq <= permrsq)
		      pval++; 

		    fprintf(ofp, "%d\t%d\t%d\t%d\t%d\t%ld\t%f\n", 
			    TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, CovCols[2]+1, p, seed, permrsq);
		  }


		  // check pval after 1000 permutations, if satisfies the condition
		  // run up to permNum
		  if ( (pval/1000)<=pseq1000 && permNum>=1000 ){

		    for (p=1001;p<=permNum;p++){
		      
		      for (i=0;i<nlabels;i++){
			N[i] = Norig[i];
			indirect[i] = i;
			invalid[i] = 0;
		      }
		      // premute the trait and calculate the permuted rsq after running rpm
		      perm_trait(num_of_rows, &tTD, TraitTest[t], &trait, N, nlabels); 
		      
		      permrsq = run_rpm(trait, N, nlabels, indirect, invalid, alpha, single, &sum, &sumsq);
		      
		      if (origrsq <= permrsq)
			pval++; 
		      
		      fprintf(ofp, "%d\t%d\t%d\t%d\t%d\t%ld\t%f\n", 
			      TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, CovCols[2]+1, p, seed, permrsq);
		    }
		  }
		}
	      }

	      fprintf(ofp_table, "%d\t%d\t%d\t%d", 
		      TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, CovCols[2]+1);
	      fprintf(ofp_table, "\t%f", origrsq);
	      // add the p value to table1 output
	      if (origrsq == 0)
		fprintf(ofp_table, "\tPval=1\n");
	      else {
		if (pval == 0)
		  fprintf(ofp_table, "\tPval=<%f\n", 1/((float) (p-1)));
		else
		  fprintf(ofp_table, "\tPval=%f\n",pval/(p-1));
	      }
	  
	      // add the grouping assignments to table1 output
	      for (i=0;i<MAXLABELS && outtable[i].group != 0;i++)
		fprintf(ofp_table, "%d %s %d\n",outtable[i].group, outtable[i].label, outtable[i].num);
	      fprintf(ofp_table, "\n");
	    }
	  }
	}
      }
    }
  }

  
  if (ncomb == 2){ // when looking at two loci at a time

    // the block of code creates the labels and then runs rpm for each combination specified.
    ind=0;
    for (t=0;t<tCounter;t++){
      for (g=0;g<covGroups;g++){
	for (l1=0;l1<(cCounter[g]-1);l1++){
	  for (l2=(l1+1);l2<cCounter[g];l2++){
	    CovCols[0] = CovTest[g][l1]; // the two cov
	    CovCols[1] = CovTest[g][l2]; // to be tested
 
	    sprintf(s,"%s", CG[i][CovCols[0]].str);
	
	    make_label(&trait, &N, &label, &nlabels, &indirect, &invalid, num_of_rows, TD, TraitTest[t], 
		       CG, CovCols, 2, sMissValue);
  
	    for (i=0;i<nlabels;i++)
	      Norig[i] = N[i]; // an array with the original information 
	    
	    origrsq = run_rpm(trait, N, nlabels, indirect, invalid, alpha, single, &sum, &sumsq);
	    
	    
	    get_output(&outtable, indirect, label, nlabels, Norig);
	    
	    
	    for (i=0;i<num_of_rows;i++){
	      for (n=0;n<nlabels;n++){
		sprintf(s,"%s", CG[i][CovCols[0]].str); 
		sprintf(s, "%s_%s", s, CG[i][CovCols[1]].str);
		
		if (strcmp(outtable[n].label,s) == 0){
		  outtable2[i][ind] = outtable[n].group;
		}
	      }
	    }
	    ind++;
	    
	    // print to file the output and table1
	    fprintf(ofp, "%d\t%d\t%d\t0\t%ld\t%f\n", TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, 
		    seed, origrsq);

	    fprintf(ofp_table, "T=%d\tL1=%d\tL2=%d\n", TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1);
	    numMissing = num_of_rows;
	    for (i=0, n=0;i<nlabels;i++) {
	      if (invalid[i] == 0){
		n++;
		numMissing-=N[i];
		var = (sumsq[i]-(pow(sum[i],2)/N[i]))/(N[i]-1);
		if (single==0 && N[i]==1)
		  fprintf(ofp_table, "G%d: mean=%f  sd=%f  n=%d  (ignored in rsq)\n",
			  n,sum[i]/N[i],sqrt(var),N[i]);
		else 
		  fprintf(ofp_table, "G%d: mean=%f  sd=%f  n=%d  (included in rsq)\n",
			  n,sum[i]/N[i],sqrt(var),N[i]);
	      }
	    }
	    fprintf(ofp_table, "Missing=%d\n",numMissing);

	    pval = 0;
	    if (origrsq != 0 || zerorsq == 0){

	      // run 100 permutations
	      for (p=1;p<=permNum && p<=100;p++){
		
		for (i=0;i<nlabels;i++){
		  N[i] = Norig[i];
		  indirect[i] = i;
		  invalid[i] = 0;
		}
		// premute the trait and calculate the permuted rsq after running rpm
		perm_trait(num_of_rows, &tTD, TraitTest[t], &trait, N, nlabels); 
		
		permrsq = run_rpm(trait, N, nlabels, indirect, invalid, alpha, single, &sum, &sumsq);
		
		if (origrsq <= permrsq)
		  pval++; 
		
		fprintf(ofp, "%d\t%d\t%d\t%d\t%ld\t%f\n", 
			TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, p, seed, permrsq);
		
	      }
	      
	      // check pval after 100 permutations, if satisfies the condition
	      // run up to 1000 permutations
	      if ( (pval/100)<=pseq100 && permNum>=100 ){
		
		for (p=101;p<=permNum && p<=1000;p++){
		  
		  for (i=0;i<nlabels;i++){
		    N[i] = Norig[i];
		    indirect[i] = i;
		    invalid[i] = 0;
		  }
		  // premute the trait and calculate the permuted rsq after running rpm
		  perm_trait(num_of_rows, &tTD, TraitTest[t], &trait, N, nlabels); 
		  
		  permrsq = run_rpm(trait, N, nlabels, indirect, invalid, alpha, single, &sum, &sumsq);
		  
		  if (origrsq <= permrsq)
		    pval++; 
		  
		  fprintf(ofp, "%d\t%d\t%d\t%d\t%ld\t%f\n", 
			  TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, p, seed, permrsq);
		  
		}
		

		// check pval after 1000 permutations, if satisfies the condition
		// run up to permNum permutations.
		if ( (pval/1000)<=pseq1000 && permNum>=1000 ){
		  
		  for (p=1001;p<=permNum;p++){
		    
		    for (i=0;i<nlabels;i++){
		      N[i] = Norig[i];
			indirect[i] = i;
			invalid[i] = 0;
		    }
		    // premute the trait and calculate the permuted rsq after running rpm
		    perm_trait(num_of_rows, &tTD, TraitTest[t], &trait, N, nlabels); 
		      
		    permrsq = run_rpm(trait, N, nlabels, indirect, invalid, alpha, single, &sum, &sumsq);
		    
		    if (origrsq <= permrsq)
			pval++; 
		    
		    fprintf(ofp, "%d\t%d\t%d\t%d\t%ld\t%f\n", 
			    TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1, p, seed, permrsq);
		  }
		}
	      }
	    }

	    fprintf(ofp_table, "%d\t%d\t%d", TraitTest[t]+1, CovCols[0]+1, CovCols[1]+1);
	    // add the p value to table1 output
	      fprintf(ofp_table, "\t%f", origrsq);
	    if (origrsq == 0)
	      fprintf(ofp_table, "\tPval=1\n");
	    else {
	      if (pval == 0)
		fprintf(ofp_table, "\tPval=<%f\n", 1/((float)(p-1)));
	      else
		fprintf(ofp_table, "\tPval=%f\n",pval/(p-1));
	    }
	    
	    // add the grouping assignments to table1 output
	    for (i=0;i<MAXLABELS && outtable[i].group != 0;i++)
	      fprintf(ofp_table, "%d %s %d\n",outtable[i].group, outtable[i].label, outtable[i].num);
	    fprintf(ofp_table, "\n");
	  }
	}
      }
    }
  }
  

  // write table2 to a file
  fprintf(ofp_table2,"ID "); 
  for (t=0;t<tNum;t++)
      fprintf(ofp_table2, "T%d ",t+1);   
  for (l1=0;l1<(cNum+gNum);l1++)
      fprintf(ofp_table2, "L%d ",l1+1);

  if (ncomb == 2){
    for (t=0;t<tCounter;t++)
      for (g=0;g<covGroups;g++)
	for (l1=0;l1<(cCounter[g]-1);l1++)
	  for (l2=(l1+1);l2<cCounter[g];l2++)
	    fprintf(ofp_table2, "T%d_L%d_L%d ",TraitTest[t]+1,CovTest[g][l1]+1,CovTest[g][l2]+1);
    fprintf(ofp_table2,"\n");
  }

  if (ncomb == 3){
    for (t=0;t<tCounter;t++)
      for (g=0;g<covGroups;g++)
	for (l1=0;l1<(cCounter[g]-2);l1++)
	  for (l2=(l1+1);l2<cCounter[g]-1;l2++)
	    for (l3=(l2+1);l3<cCounter[g];l3++)
	      fprintf(ofp_table2, "T%d_L%d_L%d_L%d ",
		      TraitTest[t]+1,CovTest[g][l1]+1,CovTest[g][l2]+1,CovTest[g][l3]+1);
    fprintf(ofp_table2,"\n");
  }


  for (i=0;i<num_of_rows;i++){
    ind = 0;
    fprintf(ofp_table2, "%d\t", ID[i]);
    for (t=0;t<tNum;t++){ 
      fprintf(ofp_table2, "%lf\t",TD[i][t]);
    }
    for (l1=0;l1<(cNum+gNum);l1++)
      fprintf(ofp_table2, "%s ",CG[i][l1]);
    

    if (ncomb == 2){
      for (t=0;t<tCounter;t++){ 
	for (g=0;g<covGroups;g++){
	  for (l1=0;l1<(cCounter[g]-1);l1++){
	    for (l2=(l1+1);l2<cCounter[g];l2++){
	      
	      fprintf(ofp_table2, "%d ", outtable2[i][ind++]);
	    }
	  }
	}
      }
      fprintf(ofp_table2, "\n");
    }

    if (ncomb == 3){
      for (t=0;t<tCounter;t++){
	for (g=0;g<covGroups;g++){
	  for (l1=0;l1<(cCounter[g]-2);l1++){
	    for (l2=(l1+1);l2<cCounter[g]-1;l2++){
	      for (l3=(l2+1);l3<cCounter[g];l3++){
		
		fprintf(ofp_table2, "%d ", outtable2[i][ind++]);
	      }
	    }
	  }
	}
      }
      fprintf(ofp_table2, "\n"); 
    }
  }
	   

  fclose(ofp);
  for (i=0;i<num_of_rows;i++){
    free(TD[i]);
    free(CG[i]);
  }
  free(ID);
  free(TD);
  free(CG);
  
  return 0;
}

/**************************************************************************************
 method: read_data                                                                       
 desc:   This method reads the param and the data files and stores all the information
         in the corresponding arrays
 param:  pfile       - name of the param file
         ifile       - name of the input file
         num_of_rows - pointer that will hold the number of samples in the data
	 tNum        - pointer that will hold the number of traits in the sample
	 traitType   - pointer that will hold the trait type
	 cNum        - pointer for covariates number
	 gNum        - pointer for genotype number
	 MissValue   - represents a missing value in the data
	 permNum     - number of permutations
	 seed        - seed value read from the file
	 zerorsq     - determines if zero rsq should be skipped
	 pseq100     - sequential testing required pval after 100 permutations
	 pseq1000    - sequential testing required pval after 1000 permutations
	 alpha       - alpha value for GH algorithm based on tuckey HSD
	 ncomb       - number of loci to look at in a simulation
	 single      - remove single sample from the calculation of Rsq
	 covGroups   - determines the number of covariate groups
	 TraitTest   - determines which traits should be used
	 tCounter    - number of traits to be analyzed
	 CovTest     - determines which covariates should be used
	 cCounter    - number of covariates in each group to be analyzed
	 ID          - array of the ID numbers
	 TD          - array of the trait values
	 CG          - array of the covariates and genotype values
 return: none
 author: Tsvika Klein                                                               
 date:   06/24/04                                                                  
 **************************************************************************************/

void read_data(char pfile[], char ifile[], int *num_of_rows, int *tNum, int *traitType, int *cNum, 
               int *gNum, int *MissValue, int *permNum, long int *seed, int *zerorsq, float *pseq100,
	       float *pseq1000, double *alpha, int *ncomb, int *single, int *covGroups, int **TraitTest, 
	       int *tCounter, int **CovTest[], int **cCounter, int **ID, double **TD[], space **CG[]){


    char *line; // a temporary space to store arrays of characters
  
    FILE *ifp,  // points pointing to the input and output files
         *pfp;  // parameter file point 
       
    fpos_t pos;  // file pointer location

    /* counters for loops */
    int k=0, i=0, t=0, c=0, g=0, j=0;

    int fnum, lnum;
   
    /* the temporary places used to do string copy and concatenate */ 
    //char buffer[20],  // store number strings
    //char str[20]; // store the concatenated string

    char *dlm = " *\n"; // delimiter
    char *dlm2 = " =\n"; // delimiter
    char *p; // used for the strtok to get the traits and cov
    char *tline; // temp line
    char *cmd;
    int cmdval;


    // open parameter file for input
    pfp = fopen(pfile,"r");
   
    // allocate memory for a line
    line  = (char *) malloc(LINESIZE * sizeof(char));
    tline = (char *) malloc(LINESIZE * sizeof(char));
    
    // test if the allocation is successful
    if (line == NULL || tline == NULL)
    {
      printf("Unable to allocate space for line\n");
    }
    
    // reads parameter file 
    while (fgets(line, LINESIZE, pfp)){
      if (strcmp(line,"\n")==0) break;

      cmd = strtok(line, dlm2);

      if (strcmp(cmd,"MISS")==0)
	(*MissValue) = atoi(strtok(NULL, dlm2));

      if (strcmp(cmd,"PERM")==0)
	(*permNum) = atoi(strtok(NULL, dlm2));

      if (strcmp(cmd,"GENO")==0)
	(*gNum) = atoi(strtok(NULL, dlm2));

      if (strcmp(cmd,"COV")==0)
	(*cNum) = atoi(strtok(NULL, dlm2));

      if (strcmp(cmd,"TRAITS")==0)
	(*tNum) = atoi(strtok(NULL, dlm2));

      if (strcmp(cmd,"TTYPE")==0)
	(*traitType) = atoi(strtok(NULL, dlm2));
      
      if (strcmp(cmd,"SEED")==0)
	(*seed) = atoi(strtok(NULL, dlm2));

      if (strcmp(cmd,"ZERO")==0)
	(*zerorsq) = atoi(strtok(NULL, dlm2));

      if (strcmp(cmd,"PSEQ")==0){
	(*pseq100)  = atof(strtok(NULL, dlm2));
	(*pseq1000) = atof(strtok(NULL, dlm2));
      }

      if (strcmp(cmd,"COMB")==0)
	(*ncomb) = atoi(strtok(NULL, dlm2));

      if (strcmp(cmd,"SINGLE")==0)
	(*single) = atoi(strtok(NULL, dlm2));
 
      if (strcmp(cmd,"ALPHA")==0)
	(*alpha) = 1 - strtod(strtok(NULL, dlm2),NULL);
    }

    
    fgets(line, LINESIZE, pfp);   
    strcpy(tline,line);
 
    // reads the traits to be tested
    t = atoi(strtok(line, dlm));
    if (t == -1){ // if -1 all traits should be tested
       *TraitTest = (int *) malloc(*tNum * sizeof(int));
       for (t=0;t<*tNum;t++) 
	 (*TraitTest)[t] = t;
       *tCounter = *tNum;
    }
    else {
      *tCounter=1;
      while ((p = (char*)strtok(NULL, dlm)) != NULL){
	printf ("p=%sp\n",p);
	  *tCounter+=1;
      }
      *TraitTest = (int *) malloc(*tCounter * sizeof(int));
      
      *tCounter = 0;
      (*TraitTest)[*tCounter] = ( atoi(strtok(tline,dlm)) - 1 );
      *tCounter = 1;
      while ((p = (char*)strtok(NULL, dlm)) != NULL){
	(*TraitTest)[*tCounter] = ( atoi(p) - 1 );
	*tCounter+=1;
      }
    }

    fgetpos(pfp,&pos);
    i=0;
    while (fgets(line, LINESIZE, pfp)){

      // ignores everything after blank line
      if (strcmp(line,"\n")==0) break; 

      *covGroups = ++i;
    }  
    fsetpos(pfp,&pos);
   
    *CovTest = (int **) malloc(*covGroups * sizeof(int *));
    *cCounter = (int *) malloc(*covGroups *  sizeof(int));

    fgetpos(pfp,&pos);

    // reads covariates that should be tested
    for (i=0;i<*covGroups;i++){
      (*cCounter)[i] = 0;
      fgets(line, LINESIZE, pfp);

      strcpy(tline,line);
      fnum = atoi(strtok(tline, dlm));
    
      (*cCounter)[i]+=1;

      if (strstr(line,"*") == NULL){
	while ((p = (char*)strtok(NULL, dlm)) != NULL)
	  (*cCounter)[i]+=1;
      }
      else {
	
        lnum = atoi(strtok(NULL, dlm));
	(*cCounter)[i] = lnum - fnum + 1;
      }
    }
    
    fsetpos(pfp,&pos);

    for (i=0;i<*covGroups;i++){
      (*CovTest)[i] = (int *) malloc( (*cCounter)[i] * sizeof(int));
      fgets(line, LINESIZE, pfp);
      strcpy(tline, line);

      (*CovTest)[i][0] = atoi(strtok(tline, dlm)) - 1;
      
      for (j=1;j<(*cCounter)[i];j++){
      
	if (strstr(line,"*") == NULL){
	  p = (char *)strtok(NULL, dlm);
	  (*CovTest)[i][j] = atoi(p) - 1;
	}
	else {
	  (*CovTest)[i][j] = (*CovTest)[i][0] + j;
	} 
      }
    }
    
    
    /* Close parameter file */
    fclose(pfp); 

    
    printf("MISS=%d\nTRAITS=%d\nCOV=%d\nGENO=%d\nPERM=%d\nTTYPE=%d\nSEED=%d\nZERO=%d\nCovGroups=%d
            \nPSEQ=%f %f\nALPHA=%f\nNCOMB=%d\nSINGLE=%d\n",
	    *MissValue,*tNum, *cNum,*gNum,*permNum,*traitType,*seed,*zerorsq, *covGroups, 
	   *pseq100, *pseq1000, *alpha, *ncomb, *single);

    printf("TraitTest[0]=%d\n",(*TraitTest)[0]);
    printf("CovTest[0][0]=%d CovTest[0][1]=%d\n",(*CovTest)[0][0],(*CovTest)[0][1]);
    //printf("CovTest[1][0]=%d CovTest[1][1]=%d\n",(*CovTest)[1][0],(*CovTest)[1][1]);
    //printf("CovTest[2][0]=%d CovTest[2][1]=%d\n",(*CovTest)[2][0],(*CovTest)[2][1]);
    //printf("CovTest[2][2]=%d CovTest[2][3]=%d\n",(*CovTest)[2][2],(*CovTest)[2][3]);
    //printf("CovTest[2][4]=%d CovTest[2][5]=%d\n",(*CovTest)[2][4],(*CovTest)[2][5]);
    printf("tCounter=%d\n",*tCounter);

    /* Open file for input */
    ifp = fopen(ifile,"r");

    /* count the number of lines in the data file */ 
    while (fgets(line, LINESIZE, ifp)  != '\0')
    {
      /* Total number of the lines of the input data */ 
      *num_of_rows+=1;
     
    }
  
    /* free the memory space used for line */
    free(line);
    free(tline);
   
    
    //store "Covariates" and "Geno Loci" 
    *CG = (struct space**) malloc (*num_of_rows * sizeof(struct space*) );
    for (i=0;i<*num_of_rows;i++){
      (*CG)[i] = (struct space*) malloc( (*cNum+*gNum) * sizeof(struct space) );
    }
   
    *ID = (int *) malloc(*num_of_rows * sizeof(int)); 
    *TD = (double **) malloc(*num_of_rows * sizeof(double*)); 
 
    for (i=0;i<*num_of_rows;i++)
      (*TD)[i] = (double*) malloc(*tNum * sizeof(double));

  
    /* start reading the data from beginning */
    freopen(ifile,"r", ifp); 
 
    /* create four string arrays for ID, Traits, Covariants, Geno Loci */ 
    for(i=0; i<*num_of_rows; i++){
	fscanf(ifp,"%d", &(*ID)[i]);
	
	for(t=0;t<*tNum;t++)
	  fscanf(ifp,"%lf", &(*TD)[i][t]);
	for(c=0;c<*cNum;c++)
	  fscanf(ifp,"%s", &(*CG)[i][c]);
	for(g=0;g<*gNum;g++)
	  fscanf(ifp,"%s", &(*CG)[i][c+g]); 

	t=c=g=0;
    }    
    fclose(ifp);
}




/**************************************************************************************
 method: make_label                                                                      
 desc:   This method creates the labels. 
 param:  trait      - array that holds all the trait values for each label. The index
                      corresponds to the index of the label array.
         N          - array that holds the number of samples in each label
         label      - array that holds the labels
         nlabels    - holds the number of labels
         indirect   - indicates the cells merged together. 
         invalid    - indicates cells that are not active anymore (after merge)
         nrows      - number of samples
         TD         - array that holds the traits for each sample
         Tcol       - trait column number that is used for the labels. The number is the
                      index in the TD array.
         CG         - array that holds all the covariates
         CovCols    - array that specifies which covariates are used for the labels. The 
                      numbers correspond to the CG array. 
         nCov       - number of covariates used for the analysis
	 sMissValue - a string that represents a  missing value
 return: none
 author: Tsvika Klein                                                               
 date:   06/24/04        
***************************************************************************************/

void make_label(double **trait[], int **N, char **label[], int *nlabels, int **indirect, int **invalid, 
		int nrows, double **TD, int Tcol, space **CG, int CovCols[], int nCov, char *sMissValue ){

  char s[nCov*10+nCov-1]; //temp string for a new label
  int currlabel; //represents a label index inside a loop 
  int c,i,j,loc1=-1,loc2=-1; //counters & locations to merge
  int delrow = 0;

  *nlabels = 0;
 
  // allocate memory to all the arrays. 
  *N = (int *) malloc(MAXLABELS * sizeof(int)); 
  *indirect = (int *) malloc(MAXLABELS * sizeof(int));
  *invalid = (int *) malloc(MAXLABELS * sizeof(int));

  *trait = (double **) malloc(MAXLABELS * sizeof(double*)); 
  *label = (char **) malloc(MAXLABELS * sizeof(char*)); 

  for (i=0;i<MAXLABELS;i++){
      (*trait)[i] = (double*) malloc(nrows * sizeof(double));
      (*label)[i] = (char*) malloc( (nCov*10+nCov-1) * sizeof(char));
      (*N)[i] = 0;
      (*indirect)[i] = i;
      (*invalid)[i] = 0;
  }
 
  for (i=0;i<nrows;i++){
 
    // remove rows with missing values
    for (c=0;c<nCov;c++){
      if (strcmp(CG[i][CovCols[c]].str,sMissValue) == 0) 
	delrow=1;
    }
    if (atoi(sMissValue) == TD[i][Tcol])
      delrow = 1;

    if (delrow == 0) {
  
      // create labels
      sprintf(s, "%s", CG[i][CovCols[0]].str); 
      for (c=1; c<nCov;c++)
	sprintf(s, "%s_%s", s, CG[i][CovCols[c]].str);
      
    
      currlabel = -1;
      for (j=0;j<*nlabels;j++){
	if (strcmp(s,(*label)[j]) == 0){
	  currlabel =j;
	  break;
	}
      }
    
      if (currlabel == -1){
	currlabel = *nlabels;
	strcpy ((*label)[currlabel],s);
	(*nlabels)++;
      }
      
      //put the trait val with the right label
      (*trait)[currlabel][(*N)[currlabel]] = TD[i][Tcol];
      (*N)[currlabel]++;
    } 
    else 
      delrow = 0;
  }

}




/**************************************************************************************
 method: run_rpm                                                                       
 desc:   This method runs the rpm program. 
 param:  trait     - array that holds all the trait values for each label. The index
                     corresponds to the index of the label array.
         N         - array that holds the number of samples in each label
         nlabels   - holds the number of labels
         indirect  - indicates the cells merged together. 
         invalid   - indicates cells that are not active anymore (after merge)
	 alpha     - alpha value for the GH method based on Tuckey
	 single    - remove single sample from the calculation of Rsq
	 sum       - holds the sum of the traits for each unique label
	 sumsq     - holds the sum of the traits squared for each unique label
         
         
 return: rsq
 author: Tsvika Klein                                                               
 date:   03/25/04        
***************************************************************************************/

double run_rpm(double **trait, int *N, int nlabels, 
	       int *indirect, int *invalid, double alpha, int single, double **sum, double **sumsq){
  
  int i,j,loc1=-1,loc2=-1, tnlabels; //counters & locations to merge
  
  double rsqval;

  // holds the sum and sumsq of the traits for each unique label
  *sum   = (double *) malloc(MAXLABELS * sizeof(double)); 
  *sumsq = (double *) malloc(MAXLABELS * sizeof(double)); 
 
  // calculate sum and sumsq for each label
  for (i=0;i<nlabels;i++){
    (*sum)[i] = 0;
    (*sumsq)[i] = 0;
   
    for (j=0;j<N[i];j++){
      (*sum)[indirect[i]] += trait[i][j];
      (*sumsq)[indirect[i]] += trait[i][j]*trait[i][j];
    }
  }

  tnlabels = nlabels;
  // merge cells according to the RPM method using GH procedure
  while (1){
    if (!choose_merge(nlabels, tnlabels, *sum, *sumsq, N, invalid, &loc1, &loc2, alpha, single)) break;
    indirect[loc2] = loc1;
    invalid[loc2] = 1;
    (*sum)[loc1] += (*sum)[loc2];
    (*sumsq)[loc1] += (*sumsq)[loc2];
    N[loc1] += N[loc2];
    tnlabels--;
  }

  rsqval = rsq(nlabels, *sum, *sumsq, N, invalid, single);

  return rsqval;
}


/**************************************************************************************
 method: choose_merge                                                                       
 desc:   This method chooses which two cells to merge using the the GH procedure by
         Games and Howell for Multiple Comparison. 
 param:  nlabels  - number of unique labels in the data
         tnlabels - updated number of labels after merge
	 sum      - sum of all the traits for each of the unique labels. After merging 
                    two entries (ie. 0,4), the entry in the low index will hold its 
                    sum plus the sum of the higher index.
	 sumsq    - sum of all the traits squared in each of the unique labels. After 
	            merging two entries (ie. 0,4), the entry in the low index will 
                    hold its sumsq plus the sumsq of the higher index.
	 N        - the number of traits (samples) in each label. 
	 invalid  - after merging two cells, the highest cell becomes inactive and 
	            marked as 1.
	 loc1     - location of the low index to be merged
	 loc2     - location of the high index to be merged
	 alpha    - alpha value for the GH method based on Tuckey
	 single   - remove single sample from the calculation of Rsq
 return: 0 if no more cells need to be merged, 1 otherwise.
 author: Tsvika Klein                                                               
 date:   06/24/04        
***************************************************************************************/

int choose_merge(int nlabels, int tnlabels, double  *sum, double *sumsq, int N[], 
		 int invalid[], int *loc1, int *loc2, double alpha, int single){
  double t,GH,df,q; // necessary variables for he GH algorithm
  int i,j, rval=0; // counters and return value.
  double mean1, mean2, meandiff, min=1000, var1, var2,temp; //var for GH
  
  // required in order to get the value for q using the function qtrng. 
  long int _l0; 
  double _f0,_f1,_f2;
  double prtrng(),qtrng(),q0000;

  
  for (i=0;i<nlabels;i++){
    for (j=i+1;j<nlabels;j++){
      if (N[i] > 1 && N[j] > 1 && invalid[i] == 0 && invalid[j] == 0){
	mean1 = sum[i]/N[i];
	mean2 = sum[j]/N[j];
	meandiff=mean1-mean2;
	//printf("meandiff=%f   ",meandiff);
	var1 = (sumsq[i]-(pow(sum[i],2)/N[i]))/(N[i]-1);
	var2 = (sumsq[j]-(pow(sum[j],2)/N[j]))/(N[j]-1);
	t = meandiff / sqrt(var1/N[i] + var2/N[j]);
	df = pow((var1/N[i] + var2/N[j]),2) / ( ( pow(var1/N[i],2)/(N[i]-1) ) + 
					        ( pow(var2/N[j],2)/(N[j]-1) ) ); 
	q = qtrng(ADR(_f0,alpha),ADR(_f1,df),ADR(_f2,tnlabels),ADR(_l0,0));
	
	if (fabs(meandiff)<min && fabs(t)<(q/sqrt(2))){
	  min=fabs(meandiff);
	  *loc1=i;
	  *loc2=j;
	  rval = 1;
	}
      }
    }
  }
  return rval;
}



/**************************************************************************************
 method: rsq                                                                       
 desc:   This method returns the rsq for the model, using SSwg and SStotal for the 
         model
 param:  nlabels   - the number of samples in the datanum of unique labels in the data
	 sum       - sum of all the traits for each of the unique labels. After merging 
                     two entries (ie. 0,4), the entry in the low index will hold its sum 
                     plus the sum of the higher index.
	 sumsq     - sum of all the traits squared in each of the unique labels. After 
	             merging two entries (ie. 0,4), the entry in the low index will hold 
                     its sumsq plus the sumsq of the higher index.
	 N         - the number of traits (samples) in each label. 
	 invalid   - after merging two cells, the highest cell becomes inactive and 
	             marked as 1.
	 single    - remove single sample from the calculation of Rsq 
 return: the req.
 author: Tsvika Klein                                                               
 date:   06/24/04        
***************************************************************************************/

double rsq(int nlabels, double  *sum, double *sumsq, int N[], int invalid[], int single){
  int i; //coutner

  // eq: sum(x^2) - sum(x)^2/Ntot = sstot1 - (sstot2)^2/sstot3 
  double sstot1=0, sstot2=0, sstot3=0; 

  // sswg=SS within groups,  sstot=SS total for the model
  double sstot=0, sswg=0;

  double r2; //rsq

  for (i=0;i<nlabels;i++){
    if (invalid[i] == 0 && (N[i]>1  || single==1)){ // Modified to remove labels with one sample
      sswg += ( sumsq[i] - ( pow(sum[i],2)/N[i] ) );

      sstot1 += sumsq[i];
      sstot2 += sum[i];
      sstot3 += N[i];
    }
  }
  sstot = sstot1 - pow(sstot2,2)/sstot3;

  r2 = 1 - sswg/sstot;
  return r2;
}




/**************************************************************************************
 method: get_output                                                                     
 desc:   this method sets up the output of table1
 param:  outtable  - array that holds the information of the output for table1         
         indirect  - indicates the cells merged together.
         label     - array that holds the labels
         nlabels   - holds the number of labels
         N         - array that holds the number of samples in each label
 return: none
 author: Tsvika Klein                                                               
 date:   06/24/04        
***************************************************************************************/

void get_output(table1 **outtable, int indirect[], char *label[], int nlabels, int N[]){
  
  
  int i, j, g, gsearch, update;
  int tindirect[MAXLABELS];

  *outtable = (table1*)malloc (MAXLABELS * sizeof(table1));

 
  // modifies indirect to hold the last groupings
  for (i=0;i<nlabels;i++){
    if (indirect[i] != i){
      indirect[i] = indirect[indirect[i]];
    }
  }

  // changes the group numbers to start from 1 with increments of 1
  g = 1;
  update = 0;
  for (gsearch=0;gsearch<nlabels;gsearch++){
    for (i=0;i<nlabels;i++){
      if (indirect[i] == gsearch){
	tindirect[i] = g;
	update = 1;
      }
    }
    if (update == 1){
      g++;
      update = 0;
    }
  }

  // updates the information in outtable
  j=0;
  for (gsearch=1;gsearch<nlabels+1;gsearch++){  
    for (i=0;i<nlabels;i++){
      if (tindirect[i] == gsearch){
	(*outtable)[j].group = tindirect[i]; 
	(*outtable)[j].label = label[i];
        (*outtable)[j].num = N[i];
	j++;
      }
    }
  }
  (*outtable)[j].group = 0;
}



/**************************************************************************************
 method: comp_ints                                                                    
 desc:   this method is used for the qsort function. The qsort function will sort the 
         arrays based on the perm value
 param:  a1  - array that holds the perm values and traits         
         a2  - array that holds the perm values and traits
 return: the value of a1->perm - a2->perm 
 author: Tsvika Klein                                                               
 date:   03/25/04        
***************************************************************************************/

int comp_ints(const struct array *a1, const struct array *a2)
{
  return (a1->perm - a2->perm);
}




/**************************************************************************************
 method: perm_traits                                                                 
 desc:   this method permutes the trait values
 param:  num_of_rows  - number of samples
         TD           - array that holds all the trait values
         Tcol         - the trait column that should be tested
         trait        - array that holds all the traits in each label
         N            - array that holds the number of samples in each label
         nlabels      - number of labels
 return: none 
 author: Tsvika Klein                                                               
 date:   03/25/04        
***************************************************************************************/

void perm_trait(int num_of_rows, double **TD[], int Tcol, double **trait[], int *N, int nlabels){
  int i,n,l;
  double val;

  struct array pt[num_of_rows]; 
  
  // permute the trait
  for (i=0;i<num_of_rows;i++){
    pt[i].trait = (*TD)[i][Tcol];
    pt[i].perm = random();
  }

  // sort te array based on the new permuted value
  qsort(pt,num_of_rows,sizeof(struct array),comp_ints);
  
  // assign the new order to the trait array
  for (i=0;i<num_of_rows;i++){
    (*TD)[i][Tcol] = pt[i].trait;
  }

  //put the trait val with the right label
  i=0;
  for (l=0;l<nlabels;l++){
    for (n=0;n<N[l];n++){  
      (*trait)[l][n] = (*TD)[i++][Tcol];
    }
  }
}




















/**************************************************************************************
 methods: prtrng, qtrng, qtrng0, alnorm, normp, ppnd  
 desc:   This methods calculate the value of q for the GH algorithm. This part is not
         documented. please refer to the reference website for details. The code was
	 converted from fortran to C using forc95.
 website: http://lib.stat.cmu.edu/apstat/   
          look under program 190 for the studentized range test
 date:   03/25/04         
***************************************************************************************/


double /*FUNCTION*/ prtrng(double *q0023,double *q0021,double *q0022,
			   long int *q0020)//long int
{
       
	long int _l0,q0040,q0041,q0039,q0033,q0033_,q0040_,_i,_r;
	double _d0,alnorm(),q0027,q0043,q0024,q0034,q0025,q0030,q0028,
	  q0042,q0044,q0035,q0031,q0032,q0000,q0037,q0002[30],q0026,q0029,
	  q0001[30],q0036,q0038;
	static double q0015[4];
	static double q0003=0.00003e0;
	static double q0004=0.0001e0;
	static double q0005=0.45e0;
	static double q0006=120.0e0;
	static double q0007=0.0e0;
	static double q0008=0.2e0;
	static double q0009=0.5e0;
	static double q0010=1.0e0;
	static double q0011=2.0e0;
	static double q0012=0.193064705e0;
	static double q0013=0.293525326e0;
	static double q0014=0.39894228e0;
	static long q0016=3;
	static long q0017=15;
	static long q0018=7;
	static long q0019=15;
	static int _aini = 1;
	if( _aini ){ /*Do 1 TIME INITIALIZATIONS!*/
		q0015[0]=0.318309886e0;
		q0015[1]=-0.268132716e-2;
		q0015[2]=0.347222222e-2;
		q0015[3]=0.833333333e-1;
		_aini = 0;
		}

	q0000=q0007;
	*q0020=0;

	if(*q0021 < q0010 || *q0022 < q0011)
	  *q0020=1;
	if(*q0023 <= q0007 || *q0020 == 1)
		goto L_99;
	q0024=q0005*pow(*q0022,-q0008);
	q0025=q0009*log(*q0022);
	q0026=*q0022-q0010;
	q0027=log(*q0022*q0024*q0014);
	if(*q0021 > q0006)
		goto L_20;
	q0028=q0005*pow(*q0021,-q0009);
	q0029=*q0021*q0009;
	if(*q0021 == q0010)
		q0027=q0012;
	if(*q0021 == q0011)
		q0027=q0013;
	if(!(*q0021 == q0010 || *q0021 == q0011))
		q0027=sqrt(q0029)*q0015[0]/(q0010+((q0015[1]/q0029+q0015[2])/
		  q0029+q0015[3])/q0029);
	q0027=log(q0027**q0022*q0024*q0028);
L_20:
	q0030=q0024;
	q0002[0]=-q0010;
	q0002[q0017]=-q0010;
	q0031=q0010;
	q0032=q0010;
	for(q0033=1;q0033<=q0019;q0033++){
		q0033_=q0033-1;
		q0030-=q0024;
L_21:
		q0030=-q0030;
		q0034=q0025+q0030;
		q0035=q0007;
		if(q0032 <= q0004 && q0033 > q0018)
			goto L_26;
		q0036=q0027-q0034*q0034*q0009;
		q0037=alnorm(&q0034,ADR(_l0,1));
		q0038=alnorm(ADR(_d0,q0034-*q0023),ADR(_l0,1))-q0037;
		if(q0038 > q0007)
			q0035=exp(q0036+q0026*log(q0038));
		if(*q0021 > q0006)
			goto L_26;
		q0039=-q0017;
L_22:
		q0039+=q0017;
		for(q0040=1;q0040<=q0017;q0040++){
			q0040_=q0040-1;
			q0041=q0040+q0039;
			if(q0002[q0041-1] > q0007)
				goto L_23;
			q0042=q0028*q0040;
			if(q0040 < q0017)
				q0002[q0041]=-q0010;
			q0043=exp(q0042);
			q0002[q0041-1]=*q0023*q0043;
			q0001[q0041-1]=*q0021*(q0042+q0009-q0043*q0043*q0009);
L_23:
			q0044=q0007;
			q0038=alnorm(ADR(_d0,q0034-q0002[q0041-1]),ADR(_l0,1))-
			  q0037;
			if(q0038 > q0007)
				q0044=exp(q0036+q0001[q0041-1]+q0026*log(q0038));
			q0035+=q0044;
			if(q0044 > q0003)
				goto L_24;
			if(q0041 > q0016 || q0033 > q0018)
				goto L_25;
L_24:
			;
			}
L_25:
		q0028=-q0028;
		if(q0028 < q0007)
			goto L_22;
L_26:
		q0000+=q0035;
		if((q0033 > q0018 && q0035 <= q0004) && q0031 <= q0004)
			goto L_99;
		q0032=q0031;
		q0031=q0035;
		if(q0030 > q0007)
			goto L_21;
		}
L_99:
	return(q0000);
} /*end of function*/

double /*FUNCTION*/ qtrng(double *q0017,double *q0015,double *q0016,
	  long int *q0013)
{
	long int q0022,q0014,q0022_;
	double q0025,q0023,q0024,q0019,q0020,prtrng(),q0018,q0021,qtrng0(),
	  q0000;
	static long q0012=8;
	static double q0001=0.001e0;
	static double q0002=0.75e0;
	static double q0003=0.80e0;
	static double q0004=0.90e0;
	static double q0005=0.99e0;
	static double q0006=0.995e0;
	static double q0007=1.75e0;
	static double q0008=1.0e0;
	static double q0009=2.0e0;
	static double q0010=5.0e0;
	static double q0011=1.0e-4;

	*q0013=0;
	q0014=0;
	if(*q0015 < q0008 || *q0016 < q0009)
		*q0013=1;
	if(*q0017 < q0004 || *q0017 > q0005)
		*q0013=2;
	if(*q0013 != 0)
		goto L_99;
	q0018=qtrng0(q0017,q0015,q0016,&q0014);
	if(q0014 != 0)
		goto L_99;
	q0019=prtrng(&q0018,q0015,q0016,&q0014);
	if(q0014 != 0)
		goto L_99;
	q0000=q0018;
	if(fabs(q0019-*q0017) < q0001)
		goto L_99;
	if(q0019 > *q0017)
		q0019=q0007**q0017-q0002*q0019;
	if(q0019 < *q0017)
		q0020=*q0017+(*q0017-q0019)*(q0008-*q0017)/(q0008-q0019)*q0002;
	if(q0020 < q0003)
		q0020=q0003;
	if(q0020 > q0006)
		q0020=q0006;
	q0021=qtrng0(&q0020,q0015,q0016,&q0014);
	if(q0014 != 0)
		goto L_99;
	for(q0022=2;q0022<=q0012;q0022++){
		q0022_=q0022-1;
		q0020=prtrng(&q0021,q0015,q0016,&q0014);
		if(q0014 != 0)
			goto L_99;
		q0023=q0019-*q0017;
		q0024=q0020-*q0017;
		q0000=(q0018+q0021)/q0009;
		q0025=q0024-q0023;
		if(fabs(q0025) > q0011)
			q0000=(q0024*q0018-q0023*q0021)/q0025;
		if(fabs(q0023) < fabs(q0024))
			goto L_12;
		q0018=q0021;
		q0019=q0020;
L_12:
		if(fabs(q0019-*q0017) < q0001*q0010)
			goto L_99;
		q0021=q0000;
		}
L_99:
	if(q0014 != 0)
		*q0013=9;
	return(q0000);
} /*end of function*/

double /*FUNCTION*/ qtrng0(double *q0011,double *q0012,double *q0014,
	  long int *ifault)
{
	double _d0,ppnd(),q0013,q0000,q0010;
	static double q0001=120.0e0;
	static double q0002=0.5e0;
	static double q0003=1.0e0;
	static double q0004=4.0e0;
	static double q0005=0.8843e0;
	static double q0006=0.2368e0;
	static double q0007=1.214e0;
	static double q0008=1.208e0;
	static double q0009=1.4142e0;

	q0010=ppnd(ADR(_d0,q0002+q0002**q0011),ifault);
	if(*q0012 < q0001)
		q0010+=(q0010*q0010*q0010+q0010)/ *q0012/q0004;
	q0013=q0005-q0006*q0010;
	if(*q0012 < q0001)
		q0013+=-q0007/ *q0012+q0008*q0010/ *q0012;
	q0000=q0010*(q0013*log(*q0014-q0003)+q0009);
	return(q0000);
} /*end of function*/

double /*FUNCTION*/ alnorm(double *q0029,int *q0027)
{
	int q0026;
	double q0000,q0030,q0028;
	static double q0001=0.0e0;
	static double q0002=1.0e0;
	static double q0003=0.5e0;
	static double q0024=7.0e0;
	static double q0025=18.66e0;
	static double q0004=1.28e0;
	static double q0005=0.398942280444e0;
	static double q0006=0.39990348504e0;
	static double q0007=0.398942280385e0;
	static double q0008=5.75885480458e0;
	static double q0009=2.62433121679e0;
	static double q0010=5.92885724438e0;
	static double q0011=-29.8213557807e0;
	static double q0012=48.6959930692e0;
	static double q0013=-3.8052e-8;
	static double q0014=3.98064794e-4;
	static double q0015=-0.151679116635e0;
	static double q0016=4.8385912808e0;
	static double q0017=0.742380924027e0;
	static double q0018=3.99019417011e0;
	static double q0019=1.00000615302e0;
	static double q0020=1.98615381364e0;
	static double q0021=5.29330324926e0;
	static double q0022=-15.1508972451e0;
	static double q0023=30.789933034e0;

	q0026=*q0027;
	q0028=*q0029;
	if(q0028 >= q0001)
		goto L_10;
	q0026=!q0026;
	q0028=-q0028;
L_10:
	if(q0028 <= q0024 || (q0026 && q0028 <= q0025))
		goto L_20;
	q0000=q0001;
	goto L_40;
L_20:
	q0030=q0003*q0028*q0028;
	if(q0028 > q0004)
		goto L_30;
	q0000=q0003-q0028*(q0005-q0006*q0030/(q0030+q0008+q0011/(q0030+
	  q0009+q0012/(q0030+q0010))));
	goto L_40;
L_30:
	q0000=q0007*exp(-q0030)/(q0028+q0013+q0019/(q0028+q0014+q0020/
	  (q0028+q0015+q0021/(q0028+q0016+q0022/(q0028+q0017+q0023/(q0028+
	  q0018))))));
L_40:
	if(!q0026)
		q0000=q0002-q0000;
	return(q0000);
} /*end of function*/

void /*FUNCTION*/ normp(double *z,double *q0019,double *q0020,double *q0018)
{
	double q0021,q0017;
	static double q0000=220.2068679123761e0;
	static double q0001=221.2135961699311e0;
	static double q0002=112.0792914978709e0;
	static double q0003=33.91286607838300e0;
	static double q0004=6.373962203531650e0;
	static double q0005=.7003830644436881e0;
	static double q0006=.3526249659989109e-01;
	static double q0007=440.4137358247522e0;
	static double q0008=793.8265125199484e0;
	static double q0009=637.3336333788311e0;
	static double q0010=296.5642487796737e0;
	static double q0011=86.78073220294608e0;
	static double q0012=16.06417757920695e0;
	static double q0013=1.755667163182642e0;
	static double q0014=.8838834764831844e-1;
	static double q0015=7.071e0;
	static double q0016=2.506628274631001e0;

	q0017=fabs(*z);
	if(q0017 > 37.e0){
		*q0018=0.e0;
		if(*z > 0.e0){
			*q0019=1.e0;
			*q0020=0.e0;
			}
		else{
			*q0019=0.e0;
			*q0020=1.e0;
			}
		return;
		}
	q0021=exp(-0.5e0*pow(q0017,(double)2));
	*q0018=q0021/q0016;
	if(q0017 < q0015){
		*q0019=q0021*((((((q0006*q0017+q0005)*q0017+q0004)*q0017+q0003)*
		  q0017+q0002)*q0017+q0001)*q0017+q0000)/(((((((q0014*q0017+
		  q0013)*q0017+q0012)*q0017+q0011)*q0017+q0010)*q0017+q0009)*
		  q0017+q0008)*q0017+q0007);
		}
	else{
		*q0019=*q0018/(q0017+1.e0/(q0017+2.e0/(q0017+3.e0/(q0017+4.e0/
		  (q0017+0.65e0)))));
		}
	if(*z < 0.e0){
		*q0020=1.e0-*q0019;
		}
	else{
		*q0020=*q0019;
		*q0019=1.e0-*q0020;
		}
	return;
} /*end of function*/

void /*FUNCTION*/ nprob(double *z,double *q0024,double *q0023,double *q0022)
{
	double q0021,q0020;
	static double q0000=0.5e0;
	static double q0001=0.398942280444e0;
	static double q0002=0.399903438504e0;
	static double q0003=5.75885480458e0;
	static double q0004=29.8213557808e0;
	static double q0005=2.62433121679e0;
	static double q0006=48.6959930692e0;
	static double q0007=5.92885724438e0;
	static double q0008=0.398942280385e0;
	static double q0009=3.8052e-8;
	static double q0010=1.00000615302e0;
	static double q0011=3.98064794e-4;
	static double q0012=1.98615381364e0;
	static double q0013=0.151679116635e0;
	static double q0014=5.29330324926e0;
	static double q0015=4.8385912808e0;
	static double q0016=15.1508972451e0;
	static double q0017=0.742380924027e0;
	static double q0018=30.789933034e0;
	static double q0019=3.99019417011e0;

	q0020=fabs(*z);
	if(q0020 > 12.7e0)
		goto L_20;
	q0021=q0000**z**z;
	*q0022=exp(-q0021)*q0008;
	if(q0020 > 1.28e0)
		goto L_10;
	*q0023=q0000-q0020*(q0001-q0002*q0021/(q0021+q0003-q0004/(q0021+
	  q0005+q0006/(q0021+q0007))));
	if(*z < 0.e0)
		goto L_30;
	*q0024=1.e0-*q0023;
	return;
L_10:
	*q0023=*q0022/(q0020-q0009+q0010/(q0020+q0011+q0012/(q0020-q0013+
	  q0014/(q0020+q0015-q0016/(q0020+q0017+q0018/(q0020+q0019))))));
	if(*z < 0.e0)
		goto L_30;
	*q0024=1.e0-*q0023;
	return;
L_20:
	*q0023=0.e0;
	*q0022=0.e0;
	if(*z < 0.e0)
		goto L_30;
	*q0024=1.e0;
	return;
L_30:
	*q0024=*q0023;
	*q0023=1.e0-*q0024;
	return;
} /*end of function*/

double /*FUNCTION*/ ppnd(double *q0021,long int *q0019)
{
	double q0000,q0020,q0022;
	static double q0001=0.42e0;
	static double q0002=2.50662823884e0;
	static double q0003=-18.61500062529e0;
	static double q0004=41.39119773534e0;
	static double q0005=-25.44106049637e0;
	static double q0006=-8.47351093090e0;
	static double q0007=23.08336743743e0;
	static double q0008=-21.06224101826e0;
	static double q0009=3.13082909833e0;
	static double q0010=-2.78718931138e0;
	static double q0011=-2.29796479134e0;
	static double q0012=4.85014127135e0;
	static double q0013=2.32121276858e0;
	static double q0014=3.54388924762e0;
	static double q0015=1.63706781897e0;
	static double q0016=0.e0;
	static double q0017=1.e0;
	static double q0018=0.5e0;

	*q0019=0;
	q0020=*q0021-q0018;
	if(fabs(q0020) > q0001)
		goto L_10;
	q0022=q0020*q0020;
	q0000=q0020*(((q0005*q0022+q0004)*q0022+q0003)*q0022+q0002)/((((q0009*
	  q0022+q0008)*q0022+q0007)*q0022+q0006)*q0022+q0017);
	return(q0000);
L_10:
	q0022=*q0021;
	if(q0020 > q0016)
		q0022=q0017-*q0021;
	if(q0022 <= q0016)
		goto L_20;
	q0022=sqrt(-log(q0022));
	q0000=(((q0013*q0022+q0012)*q0022+q0011)*q0022+q0010)/((q0015*
	  q0022+q0014)*q0022+q0017);
	if(q0020 < q0016)
		q0000=-q0000;
	return(q0000);
L_20:
	*q0019=1;
	q0000=q0016;
	return(q0000);
} /*end of function*/