#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <time.h>
#include <string.h>
#include <ctype.h>
#include <malloc.h>
#define seqLength 1000

void usage();


void usage()
{
	printf("usage: simpleCons alignment cutoff\n");
	exit(1);
}

int main(int argc, char **argv)
{
	int i,j,k;
	if(argc<3) usage();
	char predFile[50],alnFile[50];
	strcpy(predFile,argv[2]);
	strcpy(alnFile,argv[1]);
	int *cutoff5,*cutoff8,index5=1,index8=1;
//	cutoff5=ivector(1,seqLength);
//	cutoff8=ivector(1,seqLength);
	cutoff5 =new int [seqLength+1];
	cutoff8 =new int [seqLength+1];
	FILE *cutoffFile;
	char line[1000];
	char seps[]=" \n";
	char *token;
	if((cutoffFile=fopen(predFile,"r"))!=NULL){
		if(fgets(line,1000,cutoffFile)!=NULL){
			token=strtok(line,seps);
			while(token!=NULL){
				if(strcmp(token,"cutoff5:")!=0){
					cutoff5[index5++]=atoi(token);
				}
				token=strtok(NULL,seps);
			}
		}
		if(fgets(line,1000,cutoffFile)!=NULL){
			token=strtok(line,seps);
			while(token!=NULL){
				if(strcmp(token,"cutoff8:")!=0){
					cutoff8[index8++]=atoi(token);
				}
				token=strtok(NULL,seps);
			}
		}
	}
	int site[1000]; char aa[1000]; double value[1000];
	int index=0;
/*	ifstream in (alnFile,ios::in);
	if(!in) 
	{ 
		cout<<"cannot open file: "<<alnFile<<endl;
		exit(1);
	}
	while(!in.eof()){
		in>>site[index]>>aa[index]>>value[index];
		cout<<site[index]<<"  "<<aa[index]<<"  "<<value[index]<<endl;
		index++;
	}
	in.close();
	for(i=0;i<index;i++) {
		printf("%5d  %c  %8.3f\n",site[i],aa[i],value[i]);
	}
*/
	FILE *in;
	if((in=fopen(alnFile,"r"))==NULL) {printf("can not open file: %s",alnFile); exit(1);}
	else {
		while(fgets(line,1000,in)!=NULL){
			token=strtok(line,seps);
			site[index]=atoi(token);
			token=strtok(NULL,seps);
			aa[index]=token[0];
			token=strtok(NULL,seps);
			value[index]=atof(token);
			index++;
		}
	}
	fclose(in);
	for(i=0;i<index;i++) {
//		printf("%5d  %c  %8.3f\n",site[i],aa[i],value[i]);
	}
	int index1[seqLength];
	for(i=0;i<seqLength;i++) index1[i]=0;
	int limit=19;
	for(i=19;i<index-1;i++) {
		if(fabs(value[i]-value[i+1])<1.0e-8) limit=i+1;
		else break;
	}
//	printf("limit: %d\n",limit);
	int indexGroup=0;
	for(i=0;i<limit;i++){
		if(fabs(value[i]-value[i+1])<1.0e-8){
			index1[indexGroup]=i+1;
		}
		else {indexGroup++; index1[indexGroup]=i+1;}
	}
//	for(i=0;i<=indexGroup;i++) printf("%d index:%d\n",i,index1[i]);

	int hit5[seqLength+1],hit8[seqLength+1];
	for(i=0;i<=indexGroup;i++){
		hit5[i]=hit8[i]=0;
		for(j=0;j<=index1[i];j++){
			for(k=1;k<index5;k++) if(site[j]==cutoff5[k]) hit5[i]++;
			for(k=1;k<index8;k++) if(site[j]==cutoff8[k]) hit8[i]++;
		}
	}
	for(i=0;i<=indexGroup;i++) printf("site: %4d    hit5: %4d    hit8: %4d\n",index1[i]+1,hit5[i],hit8[i]);
	
	double A5_5,A5_8;
	double A10_5,A10_8;
	double A15_5,A15_8;
	double A20_5,A20_8;

	for(i=0;i<=indexGroup;i++) index1[i]+=1;
	//5
	if(index1[0]>=5) A5_5=5.0*hit5[0]/index1[0];
	else{
		for(i=1;i<=indexGroup;i++) {
			if(index1[i]==5) A5_5=(float) hit5[i];
			else {
				if(index1[i]>5 && index1[i-1]<5) A5_5=1.0*(hit5[i]-hit5[i-1])*(5-index1[i-1])/(index1[i]-index1[i-1])+hit5[i-1];
			}
		}
	}
	if(index1[0]>=5) A5_8=5.0*hit8[0]/index1[0];
	else{
		for(i=1;i<=indexGroup;i++) {
			if(index1[i]==5) A5_8=(float) hit8[i];
			else {
				if(index1[i]>5 && index1[i-1]<5) A5_8=1.0*(hit8[i]-hit8[i-1])*(5-index1[i-1])/(index1[i]-index1[i-1])+hit8[i-1];
			}
		}
	}
	//10
	if(index1[0]>=10) A10_5=10.0*hit5[0]/index1[0];
	else{
		for(i=1;i<=indexGroup;i++) {
			if(index1[i]==10) A10_5=(float) hit5[i];
			else {
				if(index1[i]>10 && index1[i-1]<10) A10_5=1.0*(hit5[i]-hit5[i-1])*(10-index1[i-1])/(index1[i]-index1[i-1])+hit5[i-1];
			}
		}
	}
	if(index1[0]>=10) A10_8=10.0*hit8[0]/index1[0];
	else{
		for(i=1;i<=indexGroup;i++) {
			if(index1[i]==10) A10_8=(float) hit8[i];
			else {
				if(index1[i]>10 && index1[i-1]<10) A10_8=1.0*(hit8[i]-hit8[i-1])*(10-index1[i-1])/(index1[i]-index1[i-1])+hit8[i-1];
			}
		}
	}
	//15
	if(index1[0]>=15) A15_5=15.0*hit5[0]/index1[0];
	else{
		for(i=1;i<=indexGroup;i++) {
			if(index1[i]==15) A15_5=(float) hit5[i];
			else {
				if(index1[i]>15 && index1[i-1]<15) A15_5=1.0*(hit5[i]-hit5[i-1])*(15-index1[i-1])/(index1[i]-index1[i-1])+hit5[i-1];
			}
		}
	}
	if(index1[0]>=15) A15_8=15.0*hit8[0]/index1[0];
	else{
		for(i=1;i<=indexGroup;i++) {
			if(index1[i]==15) A15_8=(float) hit8[i];
			else {
				if(index1[i]>15 && index1[i-1]<15) A15_8=1.0*(hit8[i]-hit8[i-1])*(15-index1[i-1])/(index1[i]-index1[i-1])+hit8[i-1];
			}
		}
	}
	//20
	if(index1[0]>=20) A20_5=20.0*hit5[0]/index1[0];
	else{
		for(i=1;i<=indexGroup;i++) {
			if(index1[i]==20) A20_5=(float) hit5[i];
			else {
				if(index1[i]>20 && index1[i-1]<20) A20_5=1.0*(hit5[i]-hit5[i-1])*(20-index1[i-1])/(index1[i]-index1[i-1])+hit5[i-1];
			}
		}
	}

//	for(i=0;i<=indexGroup;i++) printf("%d  hit5: %d   hit8: %d\n",i,hit5[i],hit8[i]);

	if(index1[0]>=20) A20_8=20.0*hit8[0]/index1[0];
	else{
		for(i=1;i<=indexGroup;i++) {
			if(index1[i]==20) { A20_8=(float) hit8[i]; }
			else {
				if(index1[i]>20 && index1[i-1]<20) A20_8=1.0*(hit8[i]-hit8[i-1])*(20-index1[i-1])/(index1[i]-index1[i-1])+hit8[i-1];
			}
		}
	}

	char conFile[50];
	strcpy(conFile,predFile);
	strcat(conFile,".con");
	FILE *out; 
	out=fopen(conFile,"a+");
	fprintf(out,"%c %c %c ",alnFile[strlen(alnFile)-7],alnFile[strlen(alnFile)-6],alnFile[strlen(alnFile)-5]);
	fprintf(out,"%5.2f %5.2f %5.2f ",A5_5,A5_8,A5_8-A5_5);
	fprintf(out,"%5.2f %5.2f %5.2f ",A10_5,A10_8,A10_8-A10_5);
	fprintf(out,"%5.2f %5.2f %5.2f ",A15_5,A15_8,A15_8-A15_5);
	fprintf(out,"%5.2f %5.2f %5.2f\n",A20_5,A20_8,A20_8-A20_5);
	fclose(out);

}

