/*********************************************************************/
/*                                                                   */
/*                             skewers.c                             */
/*            Program for matrix comparison test of random           */
/*            skewers folowing Cheverud et al. (1983) and            */
/*                         Cheverud (1996).                          */
/*           written by Liam J. Revell (ljrevell@wustl.edu)          */
/*                                                                   */
/*********************************************************************/

#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)

#define ERROR 0.000000001

// function declarations

FILE *fp_read(char *, int);
FILE *fp_write(char *, int);
double ***mem_3d(int, int, int);
double **mem_alloc(int, int);
void zero_array(double **, int, int);
double *mem_1d(int);
double ran_ecu(long *);
void rescale_vector(double *, int, double);
void zero_vector(double *, int);
double vector_correl(double *, double *, int);
void write_output(double **, int, FILE *);
double gauss(long *);


// end function declarations

// functions

// function returns gaussian random number

double gauss(long *seed)
{
	double unif[2], rsq, factor;
	int i;

	while(1)
	{
		for(i=0; i<2; i++)
			unif[i]=2.0*ran_ecu(seed)-1.0;
		rsq=unif[0]*unif[0]+unif[1]*unif[1];
		if(rsq<1) break;
	}

	factor=sqrt(-2.0*log(rsq)/rsq);
	return(unif[0]*factor);
}

// function returns file pointer "r"

FILE *fp_read(char *filename, int str_size){

	FILE *fp;
	
	filename[str_size-1]='\0';
	
	fp=fopen(filename,"r");
	if(!fp){
		printf("Failed to open %s.\nExiting.\n",filename);
		exit(1);
	}
	
	return fp;
}

// function returns file pointer "w"

FILE *fp_write(char *filename, int str_size){

	FILE *fp;
	
	filename[str_size-1]='\0';
	
	fp=fopen(filename,"w");
	if(!fp){
		printf("Failed to open %s.\nExiting.\n",filename);
		exit(1);
	}
	
	return fp;
}

// memory allocation for 3D double array

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

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

// memory allocation for 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;
}

// zero array

void zero_array(double **array, int width, int length){
	int i, j;
	for(i=0;i<width;i++) for(j=0;j<length;j++) array[i][j]=0.0;
}

// memory allocation for 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;
}

// uniform random number generator

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);
}

// rescale vector to arbitrary length

void rescale_vector(double *vector, int size, double scale){
	double sqlength, length;
	int i;
	
	sqlength=0.0;
  for(i=0;i<size;i++) sqlength+=vector[i]*vector[i];
 	length=sqrt(sqlength);
  for(i=0;i<size;i++) vector[i]/=(length/scale);

}

// zero vector

void zero_vector(double *vector, int length){
	int i;
	for(i=0;i<length;i++) vector[i]=0.0;
}

// returns vector-correlation coefficient

double vector_correl(double *vector1, double *vector2, int n){

	double dot12, dot11, dot22;
	double vc;
	int j;
	
	dot12=dot11=dot22=0.0;
	
  for(j=0;j<n;j++){
	  dot12+=vector1[j]*vector2[j];
	  dot11+=vector1[j]*vector1[j];
	  dot22+=vector2[j]*vector2[j];
	}
	// r=scp/(double)n;
  // r=(scp-sum1*sum2/(double)n)/sqrt((ss1-sum1*sum1/(double)n)*(ss2-sum2*sum2/(double)n));
  vc=dot12/sqrt(dot11*dot22);

  return vc;
}

// writes matrix to file

void write_output(double **matrix, int n, FILE *fp){

	int j, k;
	fprintf(fp,"       ");
	for(j=0;j<n;j++) fprintf(fp,"\tmatrix%d ",j+1);
	fprintf(fp,"\n");
	for(j=0;j<n;j++){
		fprintf(fp,"matrix%d",j+1);
		for(k=0;k<n;k++){
			fprintf(fp,"\t%lf",matrix[j][k]);
		}
		fprintf(fp,"\n");
	}
	fprintf(fp,"\n");	

}

// main function

int main(){
	char filename[40];
	FILE *input;
	int num_matrices, num_traits;
	double ***matrices;
	int i, j, k;
	int num_skew;
	double **r;
	double *skewer;
	long seed;
	double *R1, *R2;
	int l, m;
	double **p, r_rand;
	FILE *output;
	char c;
	
	// initialize seed
	seed=-time(NULL);
	
	// clear screen
	system("clear"); // to compile for Windows, comment this line
	// system("cls"); // to compile for Linux/Unix, comment this line
	
	// get input filename
	printf("Enter name of input file containing matrices : ");
	fgets(filename,sizeof(filename),stdin);
	
	// allocate for file pointer to input file
	input=fp_read(filename,strlen(filename));
	
	// get output filename
	printf("Enter name for output file (will overwrite) : ");
	fgets(filename,sizeof(filename),stdin);
	
	// allocate for file pointer to output file
	output=fp_write(filename,strlen(filename));
	
	// get number of matrices and traits
	fscanf(input,"%d %d",&num_matrices, &num_traits);
	
	// allocate for matrices
	matrices=mem_3d(num_matrices,num_traits,num_traits);

	// get data from input file
	for(i=0;i<num_matrices;i++){
		for(j=0;j<num_traits;j++){
			for(k=0;k<num_traits;k++){
				fscanf(input,"%lf",&matrices[i][j][k]);
			}
		}
	}
	
	// get number of skewers
	printf("Enter the number of skewers : ");
	scanf("%d",&num_skew);
	
	// allocate memory
	r=mem_alloc(num_matrices, num_matrices);
	skewer=mem_1d(num_traits);
	R1=mem_1d(num_traits);
	R2=mem_1d(num_traits);
	p=mem_alloc(num_matrices, num_matrices);
	
	// zero correlation matrix
	zero_array(r,num_matrices,num_matrices);
	zero_array(p,num_matrices,num_matrices);
	
	// calculate mean vector correlations from random skewers
	printf("Calculating mean vector correlations from random skewers.\n");
	fflush(stdout);
	
	for(i=0;i<num_skew;i++){
		
		// generate a random skewer
		for(j=0;j<num_traits;j++) skewer[j]=gauss(&seed);
//skewer[0]=1.0607; skewer[1]=0.3536;

		// rescale to length 1.0
		// rescale_vector(skewer,num_traits,1.0);
		
		for(j=0;j<num_matrices;j++){
			for(k=j;k<num_matrices;k++){
				
				// zero response vectors
				zero_vector(R1,num_traits);
				zero_vector(R2,num_traits);
				
				// multiply matrices by skewer
				for(l=0;l<num_traits;l++){
					for(m=0;m<num_traits;m++){
						R1[l]+=matrices[j][l][m]*skewer[m];
						R2[l]+=matrices[k][l][m]*skewer[m];
					}
				}

				// calculate the correlation between the vectors and average
				r[j][k]+=vector_correl(R1,R2,num_traits)/(double)num_skew;
				r[k][j]=r[j][k];
				
			}
		}
	}
	printf("\n");
	if(num_matrices<=4) write_output(r,num_matrices,stdout);
	fprintf(output,"Random skewers correlation matrix:\n");
	write_output(r,num_matrices,output);
	
	printf("Calculating the significance of the observed vector correlations.\n");
	fflush(stdout);
	
	for(i=0;i<num_skew;i++){
	
		// generate two random vectors
		for(j=0;j<num_traits;j++){ 
			R1[j]=gauss(&seed);
			R2[j]=gauss(&seed);
		}
		
		// calculate the correlation between the random vectors
		r_rand=vector_correl(R1,R2,num_traits);
		
		for(k=0;k<num_matrices;k++){
			for(l=0;l<num_matrices;l++){
				if(r_rand>=(r[k][l]-ERROR)) p[k][l]+=1.0/(double)num_skew;
			}
		}
	}
	printf("\n");
	if(num_matrices<=4) write_output(p,num_matrices,stdout);
	fprintf(output,"Random skewers significances matrix:\n");
	write_output(p,num_matrices,output);
	
	printf("Done. Output (also) printed to %s.\n\n",filename);
		
	printf("Press ENTER to exit.\n");
	c=fgetc(stdin); c=fgetc(stdin);
	
	system("clear"); // to compile for Windows, comment this line	
	// system("cls"); // to compile for Linux/Unix, comment this line
	
	return 0;
}