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

                            T_method.c  
         Program for the T-method of matrix comparison 
                  following Roff et al. (1999).  
           written by Liam J. Revell <ljrevell@wustl.edu>  
		Plan 9 port by Peter A. Cejchan <cej@gli.cas.cz>
                                                                 
**********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

#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 52744
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)

// function declarations

double ran_ecu(long *seed);
double **mem_alloc(int, int);
void get_popinfo(FILE *, char *, int *);
void get_charnum(FILE *, int *);
void get_data(FILE *, int, int, int *, double **);
int *mem_1d_int(int);
void calc_cov(int, int, int, int, int *, double **, double **);
void permute(int, int, double **, long *);
void rem_means(int,int,int,int *,double **);
double *mem_1d(int);
double calc_T(int vars, double **mat1, double **mat2);
double calc_T2(int vars, double **mat1, double **mat2);
double *mem_1d(int length);

// end

// begin functions

// this function calculates the summed absolute values of the differences between two matrices
  
double calc_T(int vars, double **mat1, double **mat2){
  int i,j;
  double T;  
  T=0;
  for(i=0;i<vars;i++){
  	for(j=i;j<vars;j++){
  		T+=fabs(mat1[i][j]-mat2[i][j]);
  	}
  }
  return T;
}

// this function calculates the summed square difference between two matrices

double calc_T2(int vars, double **mat1, double **mat2){
  int i,j;
  double T2;  
  T2=0;
  for(i=0;i<vars;i++){
  	for(j=i;j<vars;j++){
  		T2+=pow(mat1[i][j]-mat2[i][j],2);
  	}
  }
  return T2;
}

// this function allocates memory to a double vector

double *mem_1d(int length){
  double *result;

  result=malloc(length*sizeof(double));
  if(!result){ printf("Memory allocation failure.\n"); exit(1);}

  return result;
}

// this function removes population means from a data array (limit, 2 populations)

void rem_means(int n1, int n2, int vars, int *id, double **data){
  int i,j;
  double *mean1, *mean2;
  mean1=mem_1d(vars);
  mean2=mem_1d(vars);  
  for(i=0;i<vars;i++) mean1[i]=mean2[i]=0.0;
  for(i=0;i<(n1+n2);i++){
    for(j=0;j<vars;j++){
      if(id[i]==1) mean1[j]+=data[i][j]/(double)n1;
      if(id[i]==2) mean2[j]+=data[i][j]/(double)n2;
    }
  }
  for(i=0;i<(n1+n2);i++){
    for(j=0;j<vars;j++){
      if(id[i]==1) data[i][j]-=mean1[j];
      if(id[i]==2) data[i][j]-=mean2[j];
    }
  }
}

// this function returns a uniform random number on the interval (0,1) using sys/time.h to seed

double ran_ecu(long *seed){
  int j;
  long k;
  static long seed2=123456789;
  static long iy=0;
  static long iv[NTAB];
  double result;
  
  if(*seed <= 0) {
      if((*seed)==0) *seed = 1;
      seed2=(*seed);
      for(j=NTAB+7; j>=0; j--){
	k=(*seed)/IQ1;
	*seed=IA1*(*seed-k*IQ1)-k*IR1;
	if(*seed<0) *seed += IM1;
	if(j<NTAB) iv[j]=*seed;
      }
      iy=iv[0];
  }

  k=(*seed)/IQ1;
  *seed=IA1*(*seed-k*IQ1)-k*IR1;
  if(*seed<0) *seed += IM1;
  k=seed2/IQ2;
  seed2=IA2*(seed2-k*IQ2)-k*IR2;
  if(seed2<0) seed2 += IM2;
  
  j=iy/NDIV;
  iy=iv[j]-seed2;
  iv[j]=*seed;
  if (iy<1) iy += IMM1;
  
  result=AM*iy;
  
  return(result);
}

// this function permutes the data array at random

void permute(int n, int vars, double **data, long *seed){
  int i,j,k;
  double *temp_arr;

  temp_arr=mem_1d(vars);

  for(i=n-1; i>0; i--){
    j=(i+1)*ran_ecu(seed);
    for(k=0;k<vars;k++) temp_arr[k]=data[i][k];
    for(k=0;k<vars;k++) data[i][k]=data[j][k];
    for(k=0;k<vars;k++) data[j][k]=temp_arr[k];
  }
  free(temp_arr);
}

// this function allocates memory to an integer vector

int *mem_1d_int(int length){
  int *result;

  result=malloc(length*sizeof(int));
  if(!result){ printf("Memory allocation failure.\n"); exit(1);}

  return result;
}

// this function gets the poupulation information from the input file

void get_popinfo(FILE *fp, char *name, int *N){
  char string[20], c;
  int i=0;
  while((string[i]=fgetc(fp))!=' ')
    i++;
  i=0;
  while((name[i]=fgetc(fp))!=' ')
    i++;
  name[i]='\0';
  i=0;
  while((string[i]=fgetc(fp))!=' ')
    i++;
  fscanf(fp,"%d",N);
}

// this function returns the number of characters in a string before a space is found

void get_charnum(FILE *fp, int *chars){
  char string[20];
  int i=0;
  while((string[i]=fgetc(fp))!=' ')
fflush(stdout);
    i++;
  fscanf(fp,"%d",chars);
}

// this function allocates memory to a double matrix

double **mem_alloc(int width, int length){
  int i;
  double **result;
  
  result=malloc(width*sizeof(double));
  if(!result){ printf("Memory allocation failure.\n"); exit(1);}

  for(i=0;i<width;i++){ 
    result[i]=malloc(length*sizeof(double));
    if(!result[i]){ printf("Memory allocation failure.\n"); exit(1);}
  }
  return result;
}

// this function gets the data from the input file

void getdata(FILE *fp, int chars, int N, int *id, double **matrix){
  int i,j;
  for(i=0;i<N;i++){
    fscanf(fp,"%d",&id[i]);
    for(j=0;j<chars;j++){
      fscanf(fp,"%lf",&matrix[i][j]);
    }
  }
}

// this function calculates a covariance matrix

void calc_cov(int popN, int totN, int pop, int chars, int *id, double **data, double **matrix){
  int i,j,k;
  double MCP, M[2];
  for(i=0;i<chars;i++){
    for(j=0;j<chars;j++){
      MCP=M[0]=M[1]=0.0;
      for(k=0;k<totN;k++){
	if(id[k]==pop){
	  MCP+=data[k][i]*data[k][j]/(double)popN;
	  M[0]+=data[k][i]/(double)popN;
	  M[1]+=data[k][j]/(double)popN;
	}
      }
      matrix[i][j]=(double)popN*(MCP-M[0]*M[1])/(double)(popN-1);
    }
  }
}

// this is the main function

int main(){
  char filename[20], pop1[10], pop2[10], yn, c;
  int i, j;
  int N1, N2, numchars, numpermute;
  int *pop;
  long int seed=-time(NULL);
  double pct=0.0, count=0.0;
  double **data, **cov_mat1, **cov_mat2;
  double P_T=0.0, P_T2=0.0;
  double T_obs, T_r;
  double T2_obs;
  FILE *fp, *dist_fp;


  printf("Enter name of file with data : ");
  fflush(stdout);
  fgets(filename, sizeof(filename), stdin);
  filename[strlen(filename)-1]='\0';

  fp=fopen(filename,"r");
  if(!fp){
    printf("Couldn't open %s\n",filename);
    exit(1);
  }
  
  printf("Record T from randomization test to output file (Y/N) : ");
  fflush(stdout);
  yn=fgetc(stdin);

  if((yn=='y')||(yn=='Y')){
  	c=fgetc(stdin);
  	printf("Enter output filename for data from randomization distribution : ");
  fflush(stdout);
  	fgets(filename,sizeof(filename),stdin);
  	filename[strlen(filename)-1]='\0';
  	
  	dist_fp=fopen(filename, "w");
  	if(!dist_fp){
  		printf("Couldn't open %s\n",filename);
  fflush(stdout);
  		exit(1);
  	}
  }

  printf("Enter desired number of randomizations for hypothesis test : ");
  fflush(stdout);
  scanf("%d",&numpermute);

  get_popinfo(fp,pop1,&N1);
  get_popinfo(fp,pop2,&N2);
  get_charnum(fp, &numchars);

  data=mem_alloc(N1+N2,numchars);
  pop=mem_1d_int(N1+N2);

  getdata(fp,numchars,N1+N2,pop,data);

  cov_mat1=mem_alloc(numchars,numchars);
  cov_mat2=mem_alloc(numchars,numchars);

  calc_cov(N1,N1+N2,1,numchars,pop,data,cov_mat1);
  calc_cov(N2,N1+N2,2,numchars,pop,data,cov_mat2);
  
  T_obs=calc_T(numchars,cov_mat1,cov_mat2);
  T2_obs=calc_T2(numchars,cov_mat1,cov_mat2);

  rem_means(N1,N2,numchars,pop,data);

  printf("Progress\n");
  printf("|0----------------------50----------------------100|\n");
  printf("|");
  fflush(stdout);

	if((yn=='y')||(yn=='Y')){
		fprintf(dist_fp,"Random.  \tT        \tT^2      \n");
		fprintf(dist_fp,"obs.     \t%lf\t%lf\n",T_obs,T2_obs);
	}

  for(i=2;i<=numpermute;i++){
    
    if((yn=='y')||(yn=='Y')){
    	fprintf(dist_fp,"%d",i);
    	for(j=0;j<(8-(int)log10(i));j++)
    		fprintf(dist_fp," ");
    	fprintf(dist_fp,"\t");
    }
    								  								 
    permute(N1+N2,numchars,data,&seed);

    calc_cov(N1,N1+N2,1,numchars,pop,data,cov_mat1);
    calc_cov(N2,N1+N2,2,numchars,pop,data,cov_mat2);

    T_r=calc_T(numchars,cov_mat1,cov_mat2);
    if(T_obs>T_r) P_T+=1.0/(double)numpermute;
    
    if((yn=='y')||(yn=='Y'))
    	fprintf(dist_fp,"%lf\t",T_r);

    T_r=calc_T2(numchars,cov_mat1,cov_mat2);
    if(T2_obs>T_r) P_T2+=1.0/(double)numpermute;
    
    if((yn=='y')||(yn=='Y'))
    	fprintf(dist_fp,"%lf\n",T_r);    

    pct=(double)i*100.0/(double)numpermute;
    if(pct>count){
      printf("o");
      fflush(stdout);
      count+=2;
    }

  }

  printf("|\nCompleted.\n");

  printf("\nResults of Analysis : \n");
  printf("Test:    \tStatistic\tone-tailed p\n");
  printf("T        \tT   = %lf\tP = %lf\n",T_obs,P_T);
  printf("T^2      \tT^2 = %lf\tP = %lf\n",T2_obs,P_T2);
  
  printf("\nPress ENTER to exit.");
  fflush(stdout);
  
	c=fgetc(stdin); c=fgetc(stdin);

  return 0;
}
