#include <fstream.h>
#include <iomanip.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include "arrays.h"
#include "maths.h"
#define MIN(A,B) ((A <= B)? A:B)
#define PROGRAM_NAME "Jack.exe"
#define BUILD_DATE "18 December 1996"
#define VERSION "0.1"

int N,G,P;
imatrix data;
int Iterations = 5000;
iVector number_of_bands;
matrix frequency_of_bands;
Vector avg_number_bands;
iVector number_in_group;

imatrix read_data(fstream& f,int& N,int& P,int& G)
{
	int i,j;
	f >> N >> P >> G;
	imatrix data(N+1,P);
	for(i=0;i<N+1;i++)
	{
		for(j=0;j<P;j++)
		{
			f >> data[i][j];
			if(f.eof())
			{
				cerr << "Problem in read_data(): can't seem to read enough data" << endl;
				cerr << "EOF occurred before entire data array was read." << endl;
				exit(37);
			}
		}
	}
	cout << "Here's the data I've just read\n";
	cout << "\n\n-----------------------------------------------\n";
	for(i=0;i<N+1;i++)
	{
		for(j=0;j<P;j++)
		{
			cout << setw(2) << data[i][j] ;
		}
		cout << endl;
	}
	cout << "\n\n-----------------------------------------------\n";
	cout << "Number of bands = " << N << "\nNumber of lanes = " << P << "\nNumber of groups = " << G << endl;
	cout << "\n\n";

//Check to make sure the main array has only 0's and 1's
	for(i=1;i<N+1;i++)
	{
		for(j=0;j<P;j++)
		{
			if( data[i][j] <0 || data[i][j] > 1)
			{
				cerr << "Problem: your data file contains values other that 0 and 1\n";
				cerr << "For example, data["<< i << "][" << j << "] = " << data[i][j] << endl;
				cerr << "Check this.\n";
				exit(7);
			}
		}
	}
/// Check to make sure labels for groups in first row of data array are okay.
	for(j=0;j<(P);j++)
	{
		if(data[0][j] <0 || data[0][j] >(G)-1)
		{
			cerr << "Problem: your data file mentions more groups than expected\n";
			cerr << "With G = " << G << " only the values 0,...," << (G-1) << " should appear in\n";
			cerr << "the first row of data. But I found data[0][" << j << "]=" << data[0][j] << endl;
			cerr << "You'd better fix it before running the program on this data\n";
			exit(7);
		}
	}
	return data;
}
void band_statistics(void)
{
  int i,j,t,g;
  double s;
/***************************************
GROUP SIZES
*/
	iVector itemp(G);
	number_in_group=itemp;
	Vector iitemp(G+1);
	avg_number_bands = iitemp;
	for(g=0;g<G;g++) number_in_group[g]=0;
	for(j=0;j<P;j++) number_in_group[data[0][j]]++;

/***************************************
COUNT THE BANDS IN EACH LANE
*/
	iVector nbtemp(P);
	number_of_bands = nbtemp;
	for(j=0;j<P;j++)
	{
		t=0;
		for(i=1;i<N+1;i++)
		{
			t += data[i][j];
		}
		number_of_bands[j]=t;
	 }

/***************************************
FREQUENCY OF BANDS
*/

	matrix  iiitemp(G+1,N);
	frequency_of_bands =iiitemp;

	for(j=0;j<G+1;j++)
	{
		for(i=0;i<N;i++) frequency_of_bands[j][i]=0;
	}
	for(i=1;i<N+1;i++)
	{
	  for(j=0;j<P;j++)
	  {
		for(g=0;g<G;g++) if(data[0][j]==g) frequency_of_bands[g][i-1] += data[i][j];
	  }
	}

	for(i=0;i<N;i++)
	{

		//cout << setw(9) << i ;
		s=0;
		for(g=0;g<G;g++) s += frequency_of_bands[g][i] ;

		frequency_of_bands[G][i]=s;
	}



	return;
}

void similarity_tests(void)
{
	int i,j,k,g,c;
	double S,t,mean,stderr,spsquare,p,df;
	Vector means(G),errors(G),degrees(G);
	matrix similarity(P,P);
	cout << "Similarity tests via a jackknife distribution of pair-wise similarity\n";
	for(j=0;j<P;j++)
	{
		for(i=0;i<=j;i++)
		{
			t=0.0;
			for(k=1;k<N+1;k++) t += data[k][i]*data[k][j];
			if(number_of_bands[i]+number_of_bands[j] ==0)
			{
				cout << "PROBLEM: total number of bands in lanes " << i << " and " << j << " is 0" << endl;
				cout << "so the similarity measure S{" << i << "," << j <<"} isn't defined.\n";
				cout << "We'll proceed, but many results that follow are wrong.\n";
				cout << "You might want to worry why you have no bands\n";
				cout << "in one or more lanes!\n";
				similarity[j][i]=0.0;

			}
			else similarity[j][i]=2.0*t/(number_of_bands[i]+number_of_bands[j]);
		}
	}
	//The jackknife says find S_i (similarity without lane i)
	// set \phi_i = n S + (n-1) s_i
	// Take S = mean( \phi_i)
	// and stderror S to be \sqrt{ \frac{\sum (\phi_i - \bar{\phi})^2}{n(n-1)}}
	//for each group
	for(g=0;g<G;g++)
	{
		//generate a list of lanes in that group
		iVector lanes(number_in_group[g]);
		Vector Si(number_in_group[g]);
		c=0;
		for(j=0;j<P;j++) if (data[0][j]==g) lanes[c++]=j;

		//for each element in this list, find the mean of all
		 for(i=0;i<number_in_group[g];i++)
		 {
			S=0.0;
			for(j=0;j<number_in_group[g];j++)
			{
				 for(k=0;k<j;k++)
				 {
					if((j !=i) && (k != i))
					{
						if (lanes[k] < lanes[j]) S += similarity[lanes[j]][lanes[k]];
						else S+=similarity[lanes[k]][lanes[j]];
					}
				 }
			}
			Si[i] = S / (  (double) (( number_in_group[g]-1)*(number_in_group[g]-2))/2.0);

		 }
		 S=0.0;
		 for(j=0;j<number_in_group[g];j++)
		 {
			for(k=0;k<j;k++)
			{
				if (lanes[k] < lanes[j]) S += similarity[lanes[j]][lanes[k]];
				else S+=similarity[lanes[k]][lanes[j]];
			}
		 }
		 S /= (  (double) (number_in_group[g] * (number_in_group[g]-1)))/2.0;

		 mean=0.0;
		 stderr=0.0;
		 cout << " Phi values: " <<endl;
		 for(i=0;i<number_in_group[g];i++)
		 {
			mean += number_in_group[g]*S - (number_in_group[g]-1)*Si[i];
			cout <<  (number_in_group[g]*S - (number_in_group[g]-1)*Si[i]) << endl;
		 }
		 mean /= number_in_group[g];
		 for(i=0;i<number_in_group[g];i++)
		 {
			stderr +=  ((number_in_group[g]*S - (number_in_group[g]-1)*Si[i])-mean)*
			((number_in_group[g]*S - (number_in_group[g]-1)*Si[i]) -mean);
		 }
		 stderr = sqrt( stderr/ (number_in_group[g]*(number_in_group[g]-1)) );
		 cout << "Group " << g << " mean similarity= " << mean << " stderr= " << stderr << " df= " << (number_in_group[g]-1) << endl;
		 means[g]=mean;
		 errors[g]=stderr;
		 degrees[g]=number_in_group[g]-1.0;


	//pairwise within group similarities that omit this lane Call this S_i
	//
	// \phi_i = nS+(n-1)S_i
	// Find mean of the \phi_i
	}
	cout << "\n\n----------------------------------------\n\n";
	/*Here we can do some t-test*/
	for(i=1;i<G;i++)
	{
		for(j=0;j<i ;j++)
		{
			cout << "Group " << i << " vs  group "<<  j << endl;
			spsquare=(degrees[i]*(degrees[i]+1)*(errors[i]*errors[i])+  degrees[j]*(degrees[j]+1)*(errors[j]*errors[j]))/(degrees[i]+degrees[j]);
			if(spsquare > 0.0)
			{
			t= (means[i]-means[j])/sqrt(spsquare*(1.0/degrees[i] + 1.0/degrees[j]));
			df= degrees[i]+degrees[j];
			p=betai(0.5*df,0.5,df/(df+t*t));
			cout << "pooled-variance   t = " << t << " df = " << df << " two-sided p = "<< p<< endl ;
			}
			else{
			if( means[i] != means[j])
			cout << "pooled-variance: Sum of squares is 0.0, so differences appear significant at p= 0.000\n";
			else cout << "pooled-variance sum of squares and mean differences are 0.0" << endl;
			}

			df = (errors[i]*errors[i] + errors[j]*errors[j]);
			df= df*df;
			if( df > 0.0)
			{
			df = df/( pow(errors[i],4)/degrees[i] + pow(errors[j],4)/degrees[j]);

			t = (means[i]-means[j])/sqrt( errors[i]*errors[i] + errors[j]*errors[j]);
			p=betai(0.5*df,0.5,df/(df+t*t));
			cout << "different-variance   t = " << t << " df = " << df << " two-sided p = "<< p<< endl ;
			}
			else{
			if( means[i] != means[j])  cout << "different-variance: Sum of squares is 0.0, so differences appear significant at p=0.000\n";
         else cout << "different-variance sum of squares and mean differences are 0.0" << endl;
			}
		}
	}


}



void process_command_line(int argc,char *argv[])
{

	cout << PROGRAM_NAME << " version " << VERSION<< ", built  "<< BUILD_DATE<<endl;
	if(argc < 2)
	{
		cout << "\n----------------------------------------------------------\n";
		cout << "Use this program by giving the command "<< endl;
		cout << argv[0] << " <datafilename>" << endl;
		cout << "\n----------------------------------------------------------\n";
		exit( 0);
	}
	fstream fin(argv[1],ios::in);

	if( fin.fail())
	{
		cout << "Sorry: can't seem to open file " << argv[1] << endl;;
		exit(11);
	}
	if(argc==3)
	{
		if(isdigit(argv[2][0])) Iterations = atoi(argv[2]);
		else
		{
			cerr << "You seem to be trying to set the number of iterations, but\n";
			cerr << "I can't treat the third argument, " << argv[2] << ", on the command line as an integer\n";
		}
	}
	data=read_data(fin,N,P,G);
	fin.close();
	return;
}
int main(int argc, char *argv[])
{
	cout << setiosflags(ios::showpoint| ios::fixed| ios::left) ;
	cout << setprecision(4);

	set_seed();

	process_command_line(argc,argv) ;



	band_statistics();

	similarity_tests()  ;
   return 1;
}