/* ---------------------------------------------------------------------------

Hypothesis test utility for use with IM surface files (*.surf).
Code written by C R Young
Last updated: 3 Nov 2005

    Copyright (C) 2005  C. R. Young

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

Input files:
"filelist.txt"			contains the names of the files that you want to read
outfile.txt_0000.surf		all of the surface files listed in "filelist.txt"
...

Useage:
IM_hypotest <burnin> 

Options:
<burnin> 								integer > 0. Ignores this many samples at the beginning of files 

Examples of usage for 500000 burnin:
./IM_hypotest 500000

Compilation (linux):
gcc IM_hypotest.c -oIM_hypotest

Also compiled on Visual C++ 2005 Express Edition Beta
   
   ---------------------------------------------------------------------------*/

#include <stdlib.h>
#include <stdio.h>
#define MAXSTR 80
#define MAXFILES 50
/* ---------------------------------------------------------------------------*/
/* Runtime errors                                                             */
/*                                                                            */
/* ---------------------------------------------------------------------------*/

/* My error message */
void MyError(error_text)
char  error_text[];
{
	fprintf(stderr,"\nRY's run-time error...\n");
	fprintf(stderr,"\n\t%s\n",error_text);
	fprintf(stderr,"\n...now exiting to system...\n\n");
	exit(1);
}


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


    FILE *fptr, *Files;
    char currentfile[MAXSTR];
	char filelist[MAXFILES][MAXSTR];
	long int count_check;
	int end, num, numfiles, file, count, diff, burnin;
	float q1, q2, qA, m1, m2, t, Nem1, Nem2;
	float q1prev, q2prev, qAprev, m1prev, m2prev, tprev;
	char check;
	int *filesize;
	int Pr_theta1, Pr_m1, Pr_Nem1;
	float Pr_Q1, Pr_M1, Pr_Q2, Pr_M2, Pr_NEM1, Pr_NEM2, qodds;
	
/*********************/
/* Begin read files  */
/*********************/

        if( argc < 2 || argc > 2) {
        	printf("usage: IM_hypotest <burnin>\n");
        	exit(1) ;
        }
	
	burnin = atoi(argv[1]);
	if (burnin < 0) {
		printf ( "Burnin doesn't make sense\n");
		MyError("usage: IM_hypotest <BURNIN: integer >= 0>\n");
	}
	
	/* Read file lookup table ("filelist.txt") and keep names in filelist */
        Files = fopen ("filelist.txt", "r");
        if (!Files)
        {
            printf ( "Error opening file!\nCould not read filelist.txt...\n");
            MyError("Required file not found!\n");
        }
        end = 0;
        numfiles = 0;
        printf("\nFiles to read:\n");
        while (end == 0) {
            if(end == 0){
                fscanf(Files,"%s", &filelist[numfiles]);
                end = feof(Files);
                if(end == 0){                     
                    printf( "%s\n", filelist[numfiles]); 
                }
                numfiles += 1;
            }
        }
        printf("\n");
        numfiles -= 1;
        if(fclose(Files) != 0) MyError("Could not close file: filelist.txt.\n");
	
	/* memory allocation trap */
        if(numfiles > MAXFILES){
            printf ( "Error, too many files! (max = %d)\n", MAXFILES);
            MyError("Recompile with MAXFILES set to a higher number.\n");
        }
	
	printf("Number of files to read = %d\n\n", numfiles);
	
	/* Initializations */
	Pr_theta1 = 0;
	Pr_m1 = 0;
	Pr_Nem1 = 0;
	count = 0;	
	count_check = 0;
	end = 0;
	
	/* allocate array that holds number of lines in each file */
	filesize = (int *) malloc( (size_t) ( (numfiles) * sizeof(int) ) );

	for (file = 0; file < numfiles; file++){
		
		sprintf(currentfile, "%s", filelist[file]);	//open file out of file list (read from "filelist.txt")
        	fptr = fopen (currentfile, "r");
        	if (!fptr)
        	{
                	printf ( "Error opening file %d!\n", file);
                	MyError("Could not read infile\n");
        	}

		printf ( "Opened file %d: %s...\n", file, currentfile);
		
		/* read the first line of the file and store this as a previous value (zero freq) */
		fscanf(fptr,"%d%f%f%f%f%f%f", &num, &q1prev, &q2prev, &qAprev, &tprev, &m1prev, &m2prev);
		

		/* keeping track of the total samples in the file */
		filesize[file] = 0;
		
		while( end == 0){ 
	
			fscanf(fptr,"%d", &num);  // frequency of the line
			fscanf(fptr,"%*c");        // get rid of tab char (fixed error 10/10/05 linux-win bug)
			
			end = feof(fptr);         // check to see if end of file has been reached
			
			if( end == 0 ) {
				
				count += num; // counts total number of MCMC iterations in all files
				count_check += 1; // counts number of lines in all files for error check
				filesize[file] += num; // counts total number of MCMC iterations in current file
				
				check = getc(fptr);			// pull out a char to check for tab
				if( check != '\t'){			// 
					ungetc(check,fptr);		// put char back in stack
					fscanf(fptr,"%f", &q1);		// pull out float for update
					getc(fptr);			// pull out extra tab
				} else {
					q1 = q1prev;			// was tab so value did not change (no extra tab here)
				}
			
				check = getc(fptr);
				if( check != '\t'){
					ungetc(check,fptr);
					fscanf(fptr,"%f", &q2);
					getc(fptr);	
				} else {
					q2 = q2prev;
				}
							
				check = getc(fptr);
				if( check != '\t'){
					ungetc(check,fptr);
					fscanf(fptr,"%f", &qA);
					getc(fptr);	
				} else {
					qA = qAprev;
				}	
						
				check = getc(fptr);
				if( check != '\t'){
					ungetc(check,fptr);
					fscanf(fptr,"%f", &t);
					getc(fptr);	
				} else {
					t = tprev;
				}
				
				check = getc(fptr);	
				if( check != '\t'){
					ungetc(check,fptr);
					fscanf(fptr,"%f", &m1);
					getc(fptr);	
				} else {
					m1 = m1prev;
				}
				
				check = getc(fptr);
				if( check != '\t'){
					ungetc(check,fptr);
					fscanf(fptr,"%f", &m2);
					getc(fptr);	
				} else {
					m2 = m2prev;
				}
				

				/* assign values read this line to prev. holders for next read */	
				q1prev = q1;
				q2prev = q2;
				qAprev = qA;
				m1prev = m1;
				m2prev = m2;
				tprev = t;

				diff = count - burnin;	// diff is difference in total count so far (including current line, 'num') and burnin
				
				if(diff < num) {
					num = diff; // if cross the burnin threshold during this line, toss out some of num not satisfying burnin
				}
				
				if(diff > 0){ // check to see if past the burnin

					if(q1 > q2) {
						Pr_theta1 += num; 
					}
				
					if(m1 > m2) {
						Pr_m1 += num; 
					}

					Nem1 = q1*m1*0.5;
					Nem2 = q2*m2*0.5;
					if(Nem1 > Nem2 ) {
						Pr_Nem1 += num;
					}
				
				}

			}
		}
		
		if(fclose(fptr) != 0) MyError("Could not close file %d\n", file);
		printf ( "Closed file %d... filesize is %d\n", file, filesize[file]);
		end = 0;
	
	}
	printf ( "\n");
	printf ( "\n");
	printf ( "Total samples in files=%d\n", count);
	printf ( "burnin=%d\n", burnin);
	printf ( "Samples included in hypothesis calculation=%d\n\n", diff);	
	printf ( "\n");
	printf ( "\n");
	printf ( "Hypothesis tests:\n");
	printf ( "\n");
	
	Pr_Q1 = (float) Pr_theta1 / (float) diff;
	Pr_M1 = (float) Pr_m1 / (float) diff;	
	Pr_NEM1 = (float) Pr_Nem1 / (float) diff;

	if (Pr_Q1 > 0.5){
		Pr_Q2 = 1.0 - Pr_Q1;
		qodds = Pr_Q1/Pr_Q2;
		printf ( "Pr(q1>q2)=%f\n", Pr_Q1);
		if(Pr_Q2 > 0.0){
			printf ( "Odds Ratio=%.2f\n\n", qodds);
		} else {
			printf ( "Odds Ratio is undefined (infinity) \n\n");
		}
	} else {
		Pr_Q2 = 1.0 - Pr_Q1;
		qodds = Pr_Q2/Pr_Q1;
		printf ( "Pr(q2>q1)=%f\n", Pr_Q2);
		if(Pr_Q1 > 0.0){
			printf ( "Odds Ratio=%.2f\n\n", qodds);
		} else {
			printf ( "Odds Ratio is undefined (infinity) \n\n");
		}
	}
	
	if (Pr_M1 > 0.5){
		Pr_M2 = 1.0 - Pr_M1;
		qodds = Pr_M1/Pr_M2;
		printf ( "Pr(m1>m2)=%f\n", Pr_M1);
		if(Pr_M2 > 0.0){
			printf ( "Odds Ratio=%.2f\n\n", qodds);
		} else {
			printf ( "Odds Ratio is undefined (infinity) \n\n");
		}
	} else {
		Pr_M2 = 1.0 - Pr_M1;
		qodds = Pr_M2/Pr_M1;
		printf ( "Pr(m2>m1)=%f\n", Pr_M2);
		if(Pr_M1 > 0.0){
			printf ( "Odds Ratio=%.2f\n\n", qodds);
		} else {
			printf ( "Odds Ratio is undefined (infinity) \n\n");
		}
	}

	if (Pr_NEM1 > 0.5){
		Pr_NEM2 = 1.0 - Pr_NEM1;
		qodds = Pr_NEM1/Pr_NEM2;
		printf ( "Pr(Nem1>Nem2)=%f\n", Pr_NEM1);
		if(Pr_NEM2 > 0.0){
			printf ( "Odds Ratio=%.2f\n\n", qodds);
		} else {
			printf ( "Odds Ratio is undefined (infinity) \n\n");
		}
	} else {
		Pr_NEM2 = 1.0 - Pr_NEM1;
		qodds = Pr_NEM2/Pr_NEM1;
		printf ( "Pr(Nem2>Nem1)=%f\n", Pr_NEM2);
		if(Pr_NEM1 > 0.0){
			printf ( "Odds Ratio=%.2f\n\n", qodds);
		} else {
			printf ( "Odds Ratio is undefined (infinity) \n\n");
		}
	}
	
}