/*********************************************************************/
/*                                                                   */
/*                          multi_mantel.c                           */
/*            Program for multiple matrix regression and             */
/*           hypothesis testing with the Mantel procedure            */
/*                            (Mantel 1967).                         */
/*           written by Liam J. Revell (ljrevell@wustl.edu)          */
/*                                                                   */
/*********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <stddef.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 SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

#define TOLERANCE 0.000000001;

// function declarations

FILE *fp_read(char *, int);
double ***mem_3d(int, int, int);
double **mem_alloc(int, int);
void get_matrices(int n, int nummatrices, double ***matrix, FILE *fp);
void get_matrix(int n, double **matrix, FILE *fp);
double *mem_1d(int);
void unfold_nx(int n, int num_matrices, double ***matrices, double **column);
void unfold(int n, double **matrix, double *column);
void calc_cov(int totN, int chars, double **data, double **matrix);
double *cov_xy(int totN, int chars, double *y, double **x);
void gaussj(double **a, int n, double **b, int m);
int *mem_int_vect(int nl,int nh);
void free_ivector(int *v,int nl,int nh);
void reg_coefs(double **V, double *c, double *beta, int N);
double **copy_array(double **, int, int);
double SS_regression(double *y, int n, double **x, int p, double *beta, double beta0);
double calc_intercept(int n, double *y, int p, double **x, double *beta);
double SS_error(double *y, int n, double **x, int p, double *beta, double beta0);
double SS_total(double *y, int n);
double F_ratio(double SS_regression, double SS_error, int n, int p);
void print_anova(FILE *fp, double R, double table[3][4], int n, int p, double **t_table);
void permute(long *seed, int n, double **matrix);
int *mem_1d_int(int length);
double ran_ecu(long *seed);
double multi_Rsq(double *y, double *beta, double *cov_yx, double **cov_xx, int n, int p);
double *se_beta(int p, int n, double *y, double Rsq, double **cov_inv);
FILE *fp_write(char *, int);
int get_matrix_form(char c);
void get_data(int form, int n, int nummatrices, double **ymatrix, double ***xmatrices, FILE *fp);
void get_matrix_ld(int n, double **matrix, FILE *fp);
void get_matrices_ld(int n, int nummatrices, double ***matrix, FILE *fp);
void get_matrix_ud(int n, double **matrix, FILE *fp);
void get_matrices_ud(int n, int nummatrices, double ***matrix, FILE *fp);
double **print_predicted(int n, int numcols, double ***x, double *beta, double beta0, FILE *fp);
double **print_residual(int n, double **y, double **predicted, FILE *fp);
double sign(double num);

// end function declarations

// functions

// returns sign (-1 or +1)

double sign(double num){

	double sign;
	
	if(num<0.0)
		sign=-1.0;
	else
		sign=1.0;
		
	return sign;

}

// prints residuals

double **print_residual(int n, double **y, double **predicted, FILE *fp){

	int i, j;

	double **residual;
	
	residual=mem_alloc(n,n);
	
	for(i=0;i<n;i++){
		for(j=i+1;j<n;j++){
			residual[i][j]=y[i][j]-predicted[i][j];
			residual[j][i]=residual[i][j];
		}
	}

	fprintf(fp,"Residual matrix:\n");
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			if(i!=j)
				fprintf(fp,"%lf\t",residual[i][j]);
			else
				fprintf(fp,"--------\t");
		}
		fprintf(fp,"\n");
	}
	fprintf(fp,"\n");
}			

// print predicted matrix to file

double **print_predicted(int n, int numcols, double ***x, double *beta, double beta0, FILE *fp){
	
	int i, j, k;
	
	double **predicted;
	
	predicted=mem_alloc(n,n);
	
	for(i=0;i<n;i++){
		for(j=i+1;j<n;j++){
			predicted[i][j]=beta0;
			for(k=1;k<=numcols;k++){
				predicted[i][j]+=x[i][j][k-1]*beta[k];
			}
			predicted[j][i]=predicted[i][j];
		}
	}
	
	fprintf(fp,"Predicted matrix:\n");
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			if(i!=j)
				fprintf(fp,"%lf\t",predicted[i][j]);
			else
				fprintf(fp,"--------\t");
		}
		fprintf(fp,"\n");
	}
	fprintf(fp,"\n");
	
	return predicted;
}

// get unfolded data

void get_unfolded(int n, int nummatrices, double **ymatrix, double ***xmatrices, FILE *fp){

		printf("Not currently ready to accept matrices in that form.\n");
		exit(1);

}

// get x matrices - upper diagonal

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

	int i, j, k;
	
	for(i=0;i<nummatrices;i++)
		for(j=0;j<n;j++)
			for(k=j+1;k<n;k++){
				fscanf(fp,"%lf",&matrix[j][k][i]);
				matrix[k][j]=matrix[j][k];
			}	
}

// get single y matrix - upper diagonal

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

	 int j, k;
	 
	 for(j=0;j<n;j++)
		for(k=j+1;k<n;k++){
			fscanf(fp,"%lf",&matrix[j][k]);
			matrix[k][j]=matrix[j][k];
		}
			
} 

// get x matrices - lower diagonal

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

	int i, j, k;
	
	for(i=0;i<nummatrices;i++)
		for(j=1;j<n;j++)
			for(k=0;k<j;k++){
				fscanf(fp,"%lf",&matrix[j][k][i]);
				matrix[k][j][i]=matrix[j][k][i];
			}
				
}

// get single y matrix - lower diagonal

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

	 int j, k;
	 
	 for(j=1;j<n;j++)
		for(k=0;k<j;k++){
			fscanf(fp,"%lf",&matrix[j][k]);
			matrix[k][j]=matrix[j][k];
		}
} 

// get data from file

void get_data(int form, int n, int nummatrices, double **ymatrix, double ***xmatrices, FILE *fp){

	if(form==0){
		
		get_matrix(n, ymatrix, fp);
			
		get_matrices(n, nummatrices, xmatrices, fp);
	
	}
	else if(form==1){
	
		get_matrix_ld(n, ymatrix, fp);
			
		get_matrices_ld(n, nummatrices, xmatrices, fp);
		
	}
	else if(form==2){
	
		get_matrix_ud(n, ymatrix, fp);
			
		get_matrices_ud(n, nummatrices, xmatrices, fp);

	}
	else if(form==3){
	
		get_unfolded(n, nummatrices, ymatrix, xmatrices, fp);
		
	}
	else {
		printf("Not currently ready to accept matrices in that form.\n");
		exit(1);
	}
}

// gets matrix form, returns integer

int get_matrix_form(char c){

	int matrix_form;

	if((c=='F')||(c=='f')){
		matrix_form=0;
		c=fgetc(stdin);
	}
	else if((c=='L')||(c=='l')){
		matrix_form=1;
		c=fgetc(stdin);
	}
	else if((c=='U')||(c=='u')){
		matrix_form=2;
		c=fgetc(stdin);
	}
	else if((c=='A')||(c=='a')){
		matrix_form=3;
		c=fgetc(stdin);
	}
	else {
		printf("No option selected. Looking for full matrices.\n");
		matrix_form=0;
	}
	
	return matrix_form;
	
}

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

// returns the SE of beta

double *se_beta(int p, int n, double *y, double Rsq, double **cov_inv){

	int i;
	double ymean=0.0, vary=0.0;
	double term1, term2;

	double *SE;

	SE=mem_1d(p+1);
	
	// calculate the SS for y
	
	for(i=0;i<n;i++)
		ymean+=y[i]/n;
	
	for(i=0;i<n;i++)
		vary+=pow(y[i]-ymean,2)/(n-1);

	// calculate the numerator for the SE
	
	term1=sqrt(vary*(1-Rsq)/((double)n-p-1));
	
	for(i=1;i<=p;i++){
		term2=sqrt(cov_inv[i][i]);
		SE[i]=term1*term2;
	}
	
	return SE;
}
		

// returns the multiple regression Rsq

double multi_Rsq(double *y, double *beta, double *cov_yx, double **cov_xx, int n, int p){

	int i, j;
	double ymean=0.0, vary=0.0;
	
	double R=0.0;
	
	// calculate the SS for y
	
	for(i=0;i<n;i++)
		ymean+=y[i]/n;
	
	for(i=0;i<n;i++)
		vary+=pow(y[i]-ymean,2)/(n-1);
	
	// calculate the value of the Rsq
	
	for(i=1;i<=p;i++)
		R+=beta[i]*sqrt(cov_xx[i][i])*cov_yx[i]/sqrt(vary*cov_xx[i][i]*vary);
	
	return R;
}
// generates random real number on 0-1

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

// malloc 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;
}

// permute rows and columns together in matrix

void permute(long *seed, int n, double **matrix){
  int i,j,k,temp;
  int *pos;
  double **refmat;
  
  refmat=copy_array(matrix,n,n);
  pos=mem_1d_int(n);
  
  for(i=0;i<n;i++) pos[i]=i;
  /* Randomize the position vector */
  for(j=n-1;j>0;j--){
    k=(j+1)*ran_ecu(seed);
    temp=pos[j];
    pos[j]=pos[k];
    pos[k]=temp;
  } 

  /* Use the randomized position vector to move the rows and columns of mat1 */
  
  for(j=0;j<n;j++){
    for(k=0;k<n;k++){
      matrix[j][k]=refmat[pos[j]][pos[k]];
    }
  }
}

// print the ANOVA table to file

void print_anova(FILE *fp, double R, double table[3][4], int n, int p, double **t_table){

	int i;

 	fprintf(fp,"Summary of Fit :\n");
  fprintf(fp,"Rsquare\t%lf\n",R);
 	fprintf(fp,"\n");
  fprintf(fp,"Analysis of variance :\n");
  fprintf(fp,"Source\tDF\tSum of Sq.\tMean Sq. \tF-stat  \tProb>=F\n");
  fprintf(fp,"Model\t%d\t%lf\t%lf\t%lf\t%lf\n",p,table[0][0],table[0][1],table[0][2],table[0][3]);
  fprintf(fp,"Error\t%d\t%lf\t%lf\n",n-p-1,table[1][0],table[1][1]);
  fprintf(fp,"Total\t%d\t%lf\n",n,table[2][0]);
  fprintf(fp,"\nParameter Estimates : \n");
  fprintf(fp,"Variable\tEstimate\tSE      \tt        \tProb>=t\n");
  fprintf(fp,"Intercept\t%lf\n",t_table[0][0]);
  for(i=1;i<=p;i++)
  	fprintf(fp,"Beta_%d  \t%lf\t%lf\t%lf\t%lf\n",i,t_table[i][0],t_table[i][1],t_table[i][2],t_table[i][3]);
	fprintf(fp,"\n");

}

// calculate the F-ratio

double F_ratio(double SS_regression, double SS_error, int n, int p){

	double MS_reg, MS_err;
	
	MS_reg=SS_regression/(double)p;
	
	MS_err=SS_error/(double)(n-p-1);
	
	return MS_reg/MS_err;
}

// calculate the total sum of squares

double SS_total(double *y, int n){

	int i, j;
	double ymean=0.0;
	
	double SS=0.0;
	
	// calculate mean y
	
	for(i=0;i<n;i++)
		ymean+=y[i]/n;
		
	for(i=0;i<n;i++)
		SS+=pow(y[i]-ymean,2);	

	return SS;
}

// calculates the error sum of squares

double SS_error(double *y, int n, double **x, int p, double *beta, double beta0){

	double *predicted_y;
	int i, j;
	
	double SS=0.0;
	
	// allocate for a vector of predicted y
	
	predicted_y=mem_1d(n);
	
	// generated predicted y

	for(i=0;i<n;i++){
		predicted_y[i]=beta0;
		for(j=0;j<p;j++){
			predicted_y[i]+=beta[j+1]*x[i][j];
		}
		SS+=pow(y[i]-predicted_y[i],2);
	}
	
	free(predicted_y);
	
	return SS;
}

// calculates the intercept in the model

double calc_intercept(int n, double *y, int p, double **x, double *beta){

	double meany=0.0, *meanx;
	int i, j;
	
	double intercept;
	
	meanx=mem_1d(p);
	
	for(i=0;i<n;i++)
		meany+=y[i]/n;
	
	for(i=0;i<p;i++){
		meanx[i]=0.0;
		for(j=0;j<n;j++){
			meanx[i]+=x[j][i]/n;
		}
	}
		
	intercept=meany;
	
	for(i=1;i<=p;i++)
		intercept-=beta[i]*meanx[i-1];

	return intercept;
}

// calculates the SS of the regression

double SS_regression(double *y, int n, double **x, int p, double *beta, double beta0){

	double *predicted_y;
	int i, j;
	double ymean=0.0;
	
	double SS=0.0;

	// allocate for a vector of predicted y
	
	predicted_y=mem_1d(n);
	
	// calculate mean y
	
	for(i=0;i<n;i++)
		ymean+=y[i]/n;
	
	// generated predicted y

	for(i=0;i<n;i++){
		predicted_y[i]=beta0;
		for(j=0;j<p;j++){
			predicted_y[i]+=beta[j+1]*x[i][j];
		}
		SS+=pow(predicted_y[i]-ymean,2);
	}
	
	free(predicted_y);
	
	return SS;
}

// copies an array

double **copy_array(double **input, int rows, int cols){

	int i, j;
	double **result;
	
	result=mem_alloc(rows, cols);
	
	for(i=0;i<rows;i++){
		for(j=0;j<cols;j++){
			result[i][j]=input[i][j];
		}
	}
			
	return result;	
}

// calculate regression coefficients

void reg_coefs(double **Vinv, double *c, double *beta, int numx){

	int i, j;

	for(i=1;i<=numx;i++){
		beta[i]=0.0;
		for(j=1;j<=numx;j++)
			beta[i]+=Vinv[i][j]*c[j];
	}
	
}

// matrix inversion by Gauss-Jordan elimination

void gaussj(double **a, int n, double **b, int m){

  int *xc, *xr, *piv;
  int i, col, row, j, k, l, ll;
  double big, dum, pivinv, temp;


  xc=mem_int_vect(1,n);
  xr=mem_int_vect(1,n);
  piv=mem_int_vect(1,n);

  for(j=1;j<=n;j++) piv[j]=0;
  for(i=1;i<=n;i++){
    big=0.0;
    for(j=1;j<=n;j++)
      if(piv[j]!=1)
				for(k=1;k<=n;k++){
	  			if(piv[k]==0){
	    			if(fabs(a[j][k])>=big){
	      			big=fabs(a[j][k]);
	      			row=j;
	      			col=k;
	    			}
	  			}	
				}
    ++(piv[col]);
    
    if(row!=col){
      for(l=1;l<=n;l++) SWAP(a[row][l],a[col][l]);
      for(l=1;l<=m;l++) SWAP(b[row][l],b[col][l]);
    }
    
    xr[i]=row;
    xc[i]=col;
    
    if(a[col][col]==0.0) { 
    	printf("Error gaussj: Singular Matrix\n");
    	exit(1);
    }
    pivinv=1.0/a[col][col];
    a[col][col]=1.0;
    for(l=1;l<=n;l++) a[col][l]*=pivinv;
    for(l=1;l<=n;l++) b[col][l]*=pivinv;
    
    for(ll=1;ll<=n;ll++)
      if(ll!=col){
				dum=a[ll][col];
				a[ll][col]=0.0;
				for(l=1;l<=n;l++) a[ll][l]-=a[col][l]*dum;
				for(l=1;l<=m;l++) a[ll][l]-=b[col][l]*dum;
      }
  }

  for(l=n;l>=1;l--){
    if(xr[l]!=xc[l]);
  }
  free_ivector(piv,1,n);
  free_ivector(xr,1,n);
  free_ivector(xc,1,n);
}

// allocate an int vector with subscript range v[nl..nh]

int *mem_int_vect(int nl,int nh){

	int *v;
	v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int))-nl;
	if (v==NULL){ 
		printf("Memory allocation error.");
		exit(1);
	}
	return v;
}

// free an int vector allocated with mem_int_vect()

void free_ivector(int *v,int nl,int nh){
	free((char*) (v+nl));
}

// this function calculates a y with xi covariance vector

double *cov_xy(int totN, int chars, double *y, double **x){
  int i, k;
  double MCP, M[2];
  
  double *cvector;
  
  cvector=mem_1d(chars+1);
  
  for(i=0;i<chars;i++){
    MCP=M[0]=M[1]=0.0;
   	for(k=0;k<totN;k++){
	  	MCP+=y[k]*x[k][i]/(double)totN;
	  	M[0]+=y[k]/(double)totN;
	  	M[1]+=x[k][i]/(double)totN;
    }
   	cvector[i+1]=(double)totN*(MCP-M[0]*M[1])/(double)(totN-1);
  }
  
  return cvector;
  
}

// this function calculates a covariance matrix

void calc_cov(int totN, int chars, 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++){
	  		MCP+=data[k][i]*data[k][j]/(double)totN;
	  		M[0]+=data[k][i]/(double)totN;
	  		M[1]+=data[k][j]/(double)totN;
      }
      matrix[i+1][j+1]=(double)totN*(MCP-M[0]*M[1])/(double)(totN-1);
    }
  }
}

// unfold y matrix

void unfold(int n, double **matrix, double *column){
  int i,j,k;
  k=0;
  for(i=0;i<n;i++){
    for(j=i+1;j<n;j++){
      column[k]=matrix[i][j];
      k++;
    }
  }
}

// unfold x matrices

void unfold_nx(int n, int num_matrices, double ***matrices, double **column){
  int i,j,k,l;
  
  for(i=0;i<num_matrices;i++){
  	l=0;
  	for(j=0;j<n;j++){
    	for(k=j+1;k<n;k++){
      	column[l][i]=matrices[j][k][i];
      	l++;
      }
    }
  }
}

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

// get x matrices

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

	int i, j, k;
	
	for(i=0;i<nummatrices;i++)
		for(j=0;j<n;j++)
			for(k=0;k<n;k++)
				fscanf(fp,"%lf",&matrix[j][k][i]);
				
}

// get single y matrix

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

	 int j, k;
	 
	 for(j=0;j<n;j++)
		for(k=0;k<n;k++)
			fscanf(fp,"%lf",&matrix[j][k]);
			
} 

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

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

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

// main function

int main(){
	
	char filename[40];
	FILE *fp;
	double **ymatrix, ***xmatrices;
	int numvars, nummatrices;
	int numperms;
	int i, j, k;
	double **xi, *y;
	double **Vmatrix;
	double *cvector;
	double **vectors;
	double **Vinv;
	double *beta;
	double SS_reg;
	double intercept;
	double SS_err;
	double SS_tot;
	double F;
	double ANOVA_table[3][4];
	long seed;
	double Rsq;
	double *SEs;
	double *t;
	double **t_table;
	double Rsq0;
	char c;
	int matrix_form;
	double **predicted, *beta_original;
	double **ymatrix_original, **residual;
	int print_resid=0;
	char outfile[40];
	double t_sign;
	
	seed=-time(NULL);
	
	system("clear"); // comment line for Windows
	// system("cls"); // comment line for Linux
	
	printf("Enter input filename : ");
	fgets(filename,sizeof(filename),stdin);
	
	fp=fp_read(filename,strlen(filename));
	
	matrix_form=0;
	
	printf("Select matrix form of input file\n(F = full; L = upper; & U = upper diagonal; A = unfolded) : ");
	c=fgetc(stdin);
	
	matrix_form=get_matrix_form(c);
	
	printf("Enter output filename : "); 
	fgets(outfile,sizeof(outfile),stdin);
	
	printf("Print predicted and residual matrices to file (Y/N) : ");
	c=fgetc(stdin);
	
	if((c=='y')||(c=='Y'))
		print_resid=1;

	// get info
	
	fscanf(fp,"%d %d", &nummatrices, &numvars);
	
	nummatrices--;
	
	// allocate for matrices
	
	ymatrix=mem_alloc(numvars, numvars);
	
	xmatrices=mem_3d(numvars, numvars, nummatrices);
	
	// get data from file
	
	get_data(matrix_form, numvars, nummatrices, ymatrix, xmatrices, fp);
	
	// get number of permutations
	
	printf("Enter desired number of permutations : ");
	scanf("%d",&numperms);
	
	// memory allocate
	
	xi=mem_alloc(numvars*(numvars-1)/2, nummatrices);
	y=mem_1d(numvars*(numvars-1)/2);
	Vmatrix=mem_alloc(nummatrices+1, nummatrices+1);
	vectors=mem_alloc(nummatrices+1,nummatrices+1);
	beta=mem_1d(nummatrices+1);
	t=mem_1d(nummatrices+1);
	t_table=mem_alloc(nummatrices+1,4);
	beta_original=mem_1d(nummatrices+1);
	
	// fill vector
  for(i=1;i<=nummatrices;i++)
    vectors[i][1]=0.0;
	
	// unfold x
	
	unfold_nx(numvars, nummatrices, xmatrices, xi);
	
	// calculate the covariance matrix
	
	calc_cov(numvars*(numvars-1)/2, nummatrices, xi, Vmatrix);
	
	// copy the covariance matrix
	
	Vinv=copy_array(Vmatrix,nummatrices+1,nummatrices+1);
	
	// invert the covariance matrix
	
	gaussj(Vinv, nummatrices, vectors, 1);
	
	// run analysis
	
	for(i=0;i<numperms;i++){
	
		// unfold y
		
		unfold(numvars, ymatrix, y);
		
		// calculate the covariance between y and the xs
		
		cvector=cov_xy(numvars*(numvars-1)/2, numvars, y, xi);
		
		// calculate the partial regression coefficients
	
		reg_coefs(Vinv, cvector, beta, nummatrices);

		// calculate the intercept of the regression
		
		intercept=calc_intercept(numvars*(numvars-1)/2, y, nummatrices, xi, beta);

		// calculate the sum of squares of the regression

		SS_reg=SS_regression(y, numvars*(numvars-1)/2, xi, nummatrices, beta, intercept);
		
		// calculate the sum of squares of the residuals
		
		SS_err=SS_error(y, numvars*(numvars-1)/2, xi, nummatrices, beta, intercept);
		
		// calculate the total sum of squares
		
		SS_tot=SS_total(y, numvars*(numvars-1)/2);
		
		// calculate F-statistic
		
		F=F_ratio(SS_reg, SS_err, numvars*(numvars-1)/2, nummatrices);
		
		// calculate Rsq
		Rsq=multi_Rsq(y, beta, cvector, Vmatrix, numvars*(numvars-1)/2, nummatrices);
		
		// calculate SEs of beta coefficients and t-statistics
		
		SEs=se_beta(nummatrices, numvars*(numvars-1)/2, y, Rsq, Vinv);
		
		for(j=1;j<=nummatrices;j++)
			t[j]=beta[j]/SEs[j];
		
		// if i==0 record to ANOVA table
		
		if(i==0){
			Rsq0=Rsq;
			
			ANOVA_table[0][0]=SS_reg;
			ANOVA_table[0][1]=SS_reg/nummatrices;
			ANOVA_table[0][2]=F;
			ANOVA_table[0][3]=0.0;
			ANOVA_table[1][0]=SS_err;
			ANOVA_table[1][1]=SS_err/(numvars*(numvars-1)/2-nummatrices-1);
			ANOVA_table[2][0]=SS_tot;

			t_table[0][0]=intercept;
			for(j=1;j<=nummatrices;j++){
				t_table[j][0]=beta[j];
				t_table[j][1]=SEs[j];
				t_table[j][2]=t[j];
				t_table[j][3]=0.0;
				beta_original[j]=beta[j];
			}
			ymatrix_original=copy_array(ymatrix,numvars,numvars);
		}
		
		// add to significance levels
		
		F+=TOLERANCE;
		if(F>=ANOVA_table[0][2])
			ANOVA_table[0][3]+=(1.0/numperms);
		
		for(j=1;j<=nummatrices;j++){
			t_sign=sign(t[j]);
			t[j]*=t_sign;
			t[j]+=TOLERANCE;
			if(t[j]>=fabs(t_table[j][2]))
				t_table[j][3]+=(1.0/numperms);
		}
		
		// permute y matrix
		
		permute(&seed, numvars, ymatrix);
		
		// free memory allocated within the loop
		
		free(cvector);
	
	}
	
	fclose(fp);
	
	fp=fp_write(outfile,strlen(outfile));
	
	system("clear"); // comment line for Windows
	// system("cls"); // comment line for Linux
	
	fprintf(fp,"Results of regression analysis on %s : \n\n",filename);
	fprintf(stdout,"Results of regression analysis on %s : \n\n",filename);
	
	print_anova(stdout, Rsq0, ANOVA_table, numvars*(numvars-1)/2, nummatrices, t_table);
	
	print_anova(fp, Rsq0, ANOVA_table, numvars*(numvars-1)/2, nummatrices, t_table);
	
	if(print_resid==1){
	
		predicted=print_predicted(numvars, nummatrices, xmatrices, beta_original, t_table[0][0], fp);
	
		residual=print_residual(numvars, ymatrix_original, predicted, fp);
	
	}
	
	printf("Output also printed to %s.\nPress ENTER to exit.\n", outfile);
	
	c=fgetc(stdin); c=fgetc(stdin);
	
	fclose(fp);
	
	return 0;
}