# include <stdio.h>
# include <math.h>


int flag_nucleotide(char x);
int flag_all_nucleotides_3(char a, char b, char c);
int flag_all_nucleotides_4(char a, char b, char c, char d);
int topology3(char x1, char x2, char x3);
int topology4(char x1, char x2, char x3, char x4);

# define LINE_LENGTH 66
# define LINE_SKIP 13
# define SEGMENT_LENGTH 10

/*
Dirk Husmeier
Written in December 2001
Revised in August 2003
 */


main()
{
	FILE *outfile ;
	int i,j,k,t,iTaxa,t_TopDef;
	int N_taxa,N,N_TopDef;
	int nSegments, topo;
	int *topo_sequence, *topoDefSites_sequence;
	char strain[12];
	char x[10];
	char** DNA_alignment;

	
	/* Read in number of taxa */
	scanf("%d %d",&N_taxa,&N);
	if (N_taxa!=3 && N_taxa!=4){
		printf("Program can only deal with 3 or 4 taxa\n");
		exit();
	} 

	/* Allocate array for DNA sequence alignment */
	DNA_alignment= (char **)malloc(N_taxa*sizeof(char *));
	for (i=0; i<N_taxa; i++){
		DNA_alignment[i]= (char *)malloc(N*sizeof(char));	
	}

	/* Allocate array for topologies */
	topo_sequence= (int *)malloc(N*sizeof(int));

	/* Allocate array for topology-defining sites */
	topoDefSites_sequence= (int *)malloc(N*sizeof(int));

	
	printf("Length of DNA sequence alignment = %d bases \n",N);
	printf("Number of taxa= %d \n",N_taxa);

	/* Compute number of segments per strain to be read in */
	if ((N/SEGMENT_LENGTH)*SEGMENT_LENGTH==N)
		nSegments=N/SEGMENT_LENGTH;
	else
		nSegments=N/SEGMENT_LENGTH+1;

	
	/* Read in data */
	for (iTaxa=0; iTaxa<N_taxa; iTaxa++){
		scanf("%s",strain);
		printf("%s \n",strain);
		t=0;  /* position in the alignment */
		for (j=0; j<nSegments; j++){  
			scanf("%s",x);
			//printf("%s \n",x);
			for (k=0; k<SEGMENT_LENGTH; k++) {
				if (t<N){
					DNA_alignment[iTaxa][t]=x[k];
					//printf("%c",DNA_alignment[iTaxa][t]);
					t++;	
				}
			}
			//printf("\n");	
		}
	}

	/* Find and classify topology-defining sites for 3 sequences*/
	if (N_taxa==3){
		t_TopDef=0;
		for (t=0; t<N; t++){
			if(flag_all_nucleotides_3(DNA_alignment[0][t],DNA_alignment[1][t],DNA_alignment[2][t])==1){
				topo= topology3(DNA_alignment[0][t],DNA_alignment[1][t],DNA_alignment[2][t]);
				if (topo > 0)  {
					topo_sequence[t_TopDef]= topo;
					topoDefSites_sequence[t_TopDef]= t;
					t_TopDef++;
				}
			}
		}
	}

	/* Find and classify topology-defining sites for 4 sequences*/
	if (N_taxa==4){
		t_TopDef=0;
		for (t=0; t<N; t++){
			if(flag_all_nucleotides_4(DNA_alignment[0][t],DNA_alignment[1][t],DNA_alignment[2][t],DNA_alignment[3][t])==1){
				topo= topology4(DNA_alignment[0][t],DNA_alignment[1][t],DNA_alignment[2][t],DNA_alignment[3][t]);
				if (topo > 0)  {
					topo_sequence[t_TopDef]= topo;
					topoDefSites_sequence[t_TopDef]= t;
					t_TopDef++;
				}
			}
		}
	}

	/* Print out results to screen */
	N_TopDef=t_TopDef;
	printf("%d \n",N_TopDef);
	for (t_TopDef=0; t_TopDef<N_TopDef; t_TopDef++){
		printf("%d ",topo_sequence[t_TopDef]);
	}
	printf("\n");

	/* Print out results to file (file for HMM software) */
	outfile = fopen("topos.dat","w");
	if (outfile==NULL) {printf("Error: Cannot write to file\n"); exit(1);}
	fprintf(outfile,"T= %d \n",N_TopDef);
	for (t_TopDef=0; t_TopDef<N_TopDef; t_TopDef++){
		fprintf(outfile,"%d ",topo_sequence[t_TopDef]);
	}
	fprintf(outfile,"\n");
	fclose(outfile);

	/* Print out topologies only to file */
	outfile = fopen("topos.in","w");
	if (outfile==NULL) {printf("Error: Cannot write to file\n"); exit(1);}
	for (t_TopDef=0; t_TopDef<N_TopDef; t_TopDef++){
		fprintf(outfile,"%d ",topo_sequence[t_TopDef]);
	}
	fprintf(outfile,"\n");
	fclose(outfile);

	/* Print out results to file */
	outfile = fopen("topoDefSites.in","w");
	if (outfile==NULL) {printf("Error: Cannot write to file\n"); exit(1);}
	for (t_TopDef=0; t_TopDef<N_TopDef; t_TopDef++){
		fprintf(outfile,"%d ",topoDefSites_sequence[t_TopDef]);
	}
	fprintf(outfile,"\n");
	fclose(outfile);
	


	/* Free memory */
	for (i=0; i<N_taxa; i++){
		free((char *) DNA_alignment[i]);	
	}
	free((char**) DNA_alignment);

	free((int*) topo_sequence);
	free((int*) topoDefSites_sequence);
}



int flag_nucleotide(char x)
/* Returns 1 if site is nucleotide, or 0 if site is a gap etc. */
{
	int out;
	if (x=='A' || x=='C' || x=='G' || x=='T')
		out=1;
	else if (x=='a' || x=='c' || x=='g' || x=='t')
		out=1;
	else
		out=0;
	return out;
}

int flag_all_nucleotides_3(char a, char b, char c)
/* Returns 1 if nucleotides everywhere, or 0 otherwise */
{
	int out;
	out= flag_nucleotide(a)*flag_nucleotide(b)*flag_nucleotide(c);
	return out;
}

int flag_all_nucleotides_4(char a, char b, char c, char d)
/* Returns 1 if nucleotides everywhere, or 0 otherwise */
{
	int out;
	out= flag_nucleotide(a)*flag_nucleotide(b)*flag_nucleotide(c)*flag_nucleotide(d);
	return out;
}

int topology3(char x1, char x2, char x3){
	int topo;
	if (x1==x2 && x1!=x3)
		topo=1;
	else if (x1==x3 && x1!=x2)
		topo=2;
	else
		topo=0;
	return topo;
}

int topology4(char x1, char x2, char x3, char x4){
	int topo;
	if (x1==x2 && x3==x4 && x1!=x3)
		topo=1;
	else if (x1==x3 && x2==x4 && x1!=x2)
		topo=2;
	else if (x1==x4 && x2==x3 && x1!=x2)
		topo=3;
	else
		topo=0;
	return topo;
}