#include <iostream>
#define boolean int

//{  This program inputs a T4klist in the file
//T4klist.Thetaktree.
//This list of trees should have data in the form
//a b c d aabb abab abba
//where (a b) (c d) is assumed to be correct.  Ties will be ignored,
//as will the values of aabb etc.

//It analyzes the list for consistency, assuming that there are relatively
//few errors in the list.
//It outputs a new list 
//T4klist.ErrorCorrect
//with some "corrections" for the sake of consistency.  The error rectification
//is done in a number of "passes", which may be specified by the user.

//A file
//ErrorCorrectT4k.comments
//tells  information
//about the process.  

//Values of inconsistency are not reliable since aabb etc are ignored
//in the error correction process.
//Values of b therefore may differ slightly from those
//obtained when ties are allowed in the lists.

//Ties are not allowed in the list under this error correction process
//and also because the tree building method does not allow
//stars in the tree.


//}





const int maxleaves = 30; // {maximum number of species allowed}
const int maxreps = 100; // {maximum number of error correcting runs plus 1}
const int labellength = 50; // {length of the labels for species}
const  int maxstringlength = 1000; //{maximal length of an identifying string}

	 
int t4first[maxleaves-3][maxleaves-3][ maxleaves-3][ maxleaves-3];
 
int t4second[maxleaves-3][maxleaves-3][ maxleaves-3][ maxleaves-3]; 
int t4third[maxleaves-3][maxleaves-3][ maxleaves-3][ maxleaves-3]; 
int t4fourth[maxleaves-3][maxleaves-3][ maxleaves-3][ maxleaves-3]; 
float t4aabb[maxleaves-3][maxleaves-3][ maxleaves-3][ maxleaves-3]; 
float t4abab[maxleaves-3][maxleaves-3][ maxleaves-3][ maxleaves-3]; 
float t4abba[maxleaves-3][maxleaves-3][ maxleaves-3][ maxleaves-3]; 
char t4code[maxleaves-3][maxleaves-3][ maxleaves-3][ maxleaves-3][3]; 
	 


   int actualleaves; // {The actual number of leaves or strings.
//        It should be less than maxleaves.}    
   
ofstream    outcomment; //{the output file, containing a description of the output tree}
ifstream    int4list; // {the file containing the trees of the list of t4}
ofstream    outt4list; // {the file containing the modified trees}
char  identification[maxstringlength]; // {identifying string for the source of data}
int  stringlength; // {actual length of the identification string}
	int	testflag;  //{flag for the kind of criterion used}
boolean	writemonitor, writeperm;  


float    aabb, abab, abba, oldaabb, oldabab, oldabba;    
int	numberreversed;
float	totaldemerits;
int	actualrepcount;
boolean	immediate;
boolean	same;
int	nonstdtype;
char	ch;
	



int	flag; // {= 0 if - counts as a separate symbol, = 1 if - is to be ignored}
 char   labels[maxleaves+1][labellength];



	
//	{******* utility variables}
	int temp,  i, j, k, l, i1, j1;
	int demerit2, demerit3, demerit4;
//		{demerit2 counts the number of demerits for (ij)(kl);
//		demerit3 counts the number of demerits for (ik)(jl);
//		demerit4 counts the number of demerits for (il)(jk).}
	int jtrue,ktrue,ltrue; // {jtrue = j+1, ktrue = k+2, etc. offset j,k,l}
	int rep;
	int repcount;  //{the number of error rectifying runs}
	int newone;
	int limit;  
	boolean done;
	

	int oldcode;
	


void finishup(void);

void finishup(void)


{	
	char ch;
    outcomment.close();
    int4list.close();
    outt4list.close();
 	cout<<"Quit the application now.\n";
 	cin.get(ch);

}




boolean check (int, int, int, int);

boolean check(int a, int b, int c,int d)
//{This function returns true if (ab)(cd) is selected by the t4code[1].
//It returns false if any two are equal.
//Note that a,b,c, d are the true labels for the taxa.}
{

	int temp;
	int first, second, third, fourth;


	if (a == b)  return false;
	else if (a == c) return false;
	else if (a == d)  return false;
	else if (b == c) return false;
	else if (b == d) return false;
	else if (c == d) return false;
	else
		{
		first = a;
		second = b;
		third = c;
		fourth = d;
		
		if (second < first)
			{
			temp = first;
			first = second;
			second = temp;	
			}

		if (third < first)
			{
			temp = first;
			first = third;
			third = temp;
			}		

		if (fourth < first)
			{
			temp = first;
			first = fourth;
			fourth = temp;
			}	
				
		if (third < second)
			{
			temp = second;
			second = third;
			third = temp;
			}		

		if (fourth < second)
			{
			temp = second;
			second = fourth;
			fourth = temp;
			}		

		if (fourth < third)
			{
			temp = third;
			third = fourth;
			fourth = temp;
			}		

//		{now in order first < second < third < fourth}
		temp = t4code[first-1][second-2][third-3][ fourth-4][1];
		if (temp == 2)
			{  //{first matches second; but first = order1 in the original list}
			if ((first == a) && (second == b)) return true;
			else if ((first == b) && (second == a)) return true;
			else if ((first == c) && (second == d)) return true;
			else if ((first == d) && (second == c)) return true;
			else return false;			
			}
		else if (temp == 3)
			{
			if ((first == a) && (third == b))  return true;
			else if ((first == b) && (third == a)) return true;
			else if ((first == c) && (third == d)) return true;
			else if ((first == d) && (third == c)) return true;
			else return false;			
			}
		else
			{
			if ((first == a) && (fourth == b)) return true;
			else if ((first == b) && (fourth == a)) return true;
			else if ((first == c) && (fourth == d)) return true;
			else if ((first == d) && (fourth == c)) return true;
			else return false;			
			
			
			}   
		
		}
	 

} // {check}


void readln(void);

void readln(void)
{
	boolean done;
	char ch;
	done = false;
	while (! done)
		{
		{
		int4list.get(ch);
		if ((ch == '\n')||(ch == '\r')) done = true;
		}
		}

}  // readln

int main()
{

 writeperm = true;
 writemonitor = true;
 numberreversed = 0;

		
 
 outcomment.open("ErrorRectifyT4k.output");
 outt4list.open("T4klist.ErrorRectify");
 
    
    
	
  
 // {read in the t4list}
  
 
 	int4list.open("T4klist.Thetaktree");         			
	int4list>>actualleaves >>flag;
   	outt4list << actualleaves << ' ' << flag  << '\n';
    		
    
             if ((actualleaves > maxleaves) || (actualleaves < 0)) 
           {
            outcomment<< "Wrong number of leaves\n";
            finishup();
            }
          

	readln();
    			
//  {read in the identifying data for the source of strings}  
  stringlength = 0;
  done = false;
  while (! done)
  	{
	int4list.get(ch);
 	if ((ch == '\n') || (ch == '\r'))
 		{
 		done = true;
 		}
  	else if (stringlength >= maxstringlength) done = true;
  	else
  		{
  		identification[stringlength]= ch;
  		stringlength = stringlength + 1;
  		}
  	}			
 
  
  for (i1 = 0; i1< stringlength; i1++)
  	{	
    	outcomment <<identification[i1];	
    	cout<<identification[i1];	
    	outt4list <<identification[i1];	
   	} 
   	outcomment << " Corrections made for consistency.\n";
   	outt4list <<" Corrections made for consistency.\n";
   	cout << '\n';	
   			
//    	  	{read in the labels for species}
 
 	for (i1 = 1; i1 <= actualleaves; i1++)
 		for (j1 = 1; j1< labellength; j1++)
 			labels[i1][j1] = ' ';
  	
  	for (i1 = 1; i1<= actualleaves; i1++)
  		{
  		int4list >> i;
  		done = false;
  		j1 = 1;
  		while (! done)
  		{
  		int4list.get(ch);
  		if (ch == '\n' ) done = true;
  		else if (j1 >= labellength)  
  			{
  			done = true;
 	   readln();			
  			}
  		else
  		{
  		labels[i1][j1]=ch;
  		j1++;}
  			
  		
  		}
  		if (i != i1) 
  			{
  			outcomment<<"****** Mismatch of species in file \n" ;
  			outcomment<<"       Species should be in numerical order.\n";
  			finishup();
  			
  			}
  		}
 	for (i1 = 1; i1<= actualleaves; i1++)
 		{
 		outt4list <<i1 << ' ';
 		for (j1 = 1; j1< labellength; j1++)
 			outt4list << labels[i1][j1];
 		outt4list << '\n';
 		}   
 
    cout<<"Program has been launched.  Reading the input file .... \n";
    			
//    {read in the actual quartets}	
    
    limit = actualleaves;
    
    
	done = false;
    rep = 1;
    				    
    while (! done)
    	{
   
     	int4list>> i >> j >> k >> l >> aabb >>abab >> abba;  
    	readln();			
    	
    	if (i == 0)	
    		{
    		done = true;
    		}
    		
    	else
    	{	

     oldcode = 2;
      	
//	{sort into order i<j<k<l}
	
	if (j < i)
		{
		temp = j;
		j = i;
		i = temp;
		oldaabb  = aabb;
		oldabab = abab;
		oldabba = abba;
		aabb = oldaabb;
		abab = oldabba;
		abba = oldabab;		
		}    	

	if (k < i)
		{
		temp = k;
		k = i;
		i = temp;
		oldaabb  = aabb;
		oldabab = abab;
		oldabba = abba;
		aabb = oldabba;
		abab = oldabab;
		abba = oldaabb;	
		if (oldcode == 2) oldcode = 4;	
		}    	
    		
 	if (l < i)
		{
		temp = l;
		l = i;
		i = temp;
		oldaabb  = aabb;
		oldabab = abab;
		oldabba = abba;
		aabb = oldabab;
		abab = oldaabb;
		abba = oldabba;	
		if (oldcode == 2) oldcode = 3;
		else if (oldcode == 3) oldcode = 2;	
		}    	

 	if (k < j)
		{
		temp = k;
		k = j;
		j = temp;
		oldaabb  = aabb;
		oldabab = abab;
		oldabba = abba;
		aabb = oldabab;
		abab = oldaabb;
		abba = oldabba;	
		if (oldcode == 2) oldcode = 3;
		else if (oldcode == 3) oldcode = 2;	
		}    	

 	if (l < j)
		{
		temp = l;
		l = j;
		j = temp;
		oldaabb  = aabb;
		oldabab = abab;
		oldabba = abba;
		aabb = oldabba;
		abab = oldabab;
		abba = oldaabb;	
		if (oldcode == 2) oldcode = 4;
		else if (oldcode == 4) oldcode = 2;	
		}    	

 	if (l < k)
		{
		temp = l;
		l = k;
		k = temp;
		oldaabb  = aabb;
		oldabab = abab;
		oldabba = abba;
		aabb = oldaabb;
		abab = oldabba;
		abba = oldabab;
		if (oldcode == 3) oldcode = 4;
		else if (oldcode == 4) oldcode = 3;		
		}    	

	if (i>j)
	outcomment<<"Error in reordering of taxa " <<i << ' '<< j << ' ' << k << ' ' << l 
<<'\n'; 

	t4aabb[i-1][j-2][k-3][l-4] = aabb;
	t4abab[i-1][j-2][k-3][l-4] = abab;
	t4abba[i-1][j-2][k-3][l-4] = abba;
   	
   	t4code[i-1][j-2][k-3][l-4][0] = oldcode;
     	
    				
    				
  } //  {else}  	

    		
       	} // {while not done}
    
    outcomment << '\n';		


 outcomment<< "There are " << actualleaves << " species. \n";
 cout << "There are " << actualleaves << " species. \n";

   for (i1 = 1; i1 <= actualleaves; i1++)
   {
   outcomment << i1 << ' ';
   for (j1=1; j1< labellength; j1++)
   outcomment  << labels[i1][j1];
  outcomment << "  \n";

   };
	outcomment << '\n';
 if (flag == 1) outcomment << "Sites containing gaps for a given quartet were ignored.\n";
 	else outcomment <<"Missing symbols were treated like an extra symbol, different 
from others.\n";
 

   for (i1 = 1; i1<= actualleaves; i1++)
   {
   
     cout << i1 << ' ';
      for (j1=1; j1< labellength; j1++)
   cout<< labels[i1][j1];
   cout << '\n';
   }
   
   

 

 cout << "Maximum number of error-correcting runs was 50.\n";
 outcomment << "Maximum number of error-correcting runs was 50.\n";
 
 immediate = false;
 repcount = 50;
 nonstdtype = 20;


   
   if (actualleaves <= 0) finishup();


testflag = 1;   

    for (i = 1; i<= limit; i++)
    for (j = i+1; j<= limit; j++)
    for (k = j+1; k<= limit; k++)
    for (l = k+1; l <= limit; l++)		
		t4code[i-1][j-2][k-3][l-4][1] = t4code[i-1][j-2][k-3][l-4][0];

same = false;


//{perform tests of consistency using code[1] version in the checks}

for (rep = 1; rep <= repcount; rep++)
if (! same)
{


	totaldemerits = 0.0;


    for (i = 1; i<= limit; i++)
    for (j = i+1; j<= limit; j++)
    for (k = j+1; k<= limit; k++)
    for (l = k+1; l <= limit; l++)		
    	
    	{ // {utilize quartet i, j, k, l}
		demerit2 = 0;
		demerit3 = 0;
		demerit4 = 0;
//		{demerit2 counts the number of demerits for (ij)(kl);
//		demerit3 counts the number of demerits for (ik)(jl);
//		demerit4 counts the number of demerits for (il)(jk).}
		
		jtrue = j;
		ktrue = k; 
		ltrue = l;

		for (newone = 1; newone <=actualleaves; newone++)
		if (newone != i)
		if (newone != jtrue)
		if (newone != ktrue)
		if (newone != ltrue)
			
			

			{
			temp = 0; //{number of 100% correct checks}
	//		{Here each possibility lists all the quartets for a given placement
	//		of newone into (ij)(kl) other than (ij)(kl) itself}
		
			if (check(i,newone,jtrue,ktrue) && check(i, newone, jtrue,ltrue)
				&& check(i, newone, ktrue,ltrue) && 
check(jtrue,newone,ktrue,ltrue))
				 temp = temp + 1;
			if (temp == 0)  if (check(jtrue,newone,i,ktrue) && 
check(jtrue,newone,i,ltrue)
				&& check(i,newone,ktrue,ltrue) && 
check(jtrue,newone,ktrue,ltrue))
				 temp = temp + 1;
			if (temp == 0)  if (check(i,jtrue,ktrue,newone) && 
check(i,jtrue,ltrue,newone)
				&& check(i,newone,ktrue,ltrue) && 
check(jtrue,newone,ktrue,ltrue))
				 temp = temp + 1;
			if (temp == 0)  if (check(i,jtrue,ktrue,newone) && 
check(i,jtrue,ltrue,newone)
				&& check(i,ltrue,ktrue,newone) && 
check(jtrue,ltrue,ktrue,newone))
				 temp = temp + 1;
			if (temp == 0)  if (check(i,jtrue,ktrue,newone) && 
check(i,jtrue,ltrue,newone)
				&& check(i,ktrue,ltrue,newone) && 
check(jtrue,ktrue,ltrue,newone))
				 temp = temp + 1;	
			if (temp == 0)  demerit2 = demerit2 + 1;	
			
			temp = 0; //{number of 100% correct checks}
	//		{Here each possibility lists all the quartets for a given placement
	//		of newone into (ik)(jl) other than (ik)(jl) itself}
			
			if (check(i,newone,ktrue,jtrue) && check(i, newone, ktrue,ltrue)
				&& check(i, newone, jtrue,ltrue) && 
check(ktrue,newone,jtrue,ltrue))
				 temp = temp + 1;
			if (temp == 0)  if (check(ktrue,newone,i,jtrue) && 
check(ktrue,newone,i,ltrue)
				&& check(i,newone,jtrue,ltrue) && 
check(ktrue,newone,jtrue,ltrue))
				 temp = temp + 1;
			if (temp == 0)  if (check(i,ktrue,jtrue,newone) && 
check(i,ktrue,ltrue,newone)
				&& check(i,newone,jtrue,ltrue) && 
check(ktrue,newone,jtrue,ltrue))
				 temp = temp + 1;
			if (temp == 0)  if (check(i,ktrue,jtrue,newone) && 
check(i,ktrue,ltrue,newone)
				&& check(i,ltrue,jtrue,newone) && 
check(ktrue,ltrue,jtrue,newone))
				 temp = temp + 1;
			if (temp == 0)  if (check(i,ktrue,jtrue,newone) && 
check(i,ktrue,ltrue,newone)
				&& check(i,jtrue,ltrue,newone) && 
check(ktrue,jtrue,ltrue,newone))
				 temp = temp + 1;	
			if (temp == 0)  demerit3 = demerit3 + 1;	
			
			temp = 0; // {number of 100% correct checks}
			if (check(i,newone,ltrue,jtrue) && check(i, newone, ltrue,ktrue)
				&& check(i, newone, jtrue,ktrue) && 
check(ltrue,newone,jtrue,ktrue))
				 temp = temp + 1;
			if (temp == 0)  if (check(ltrue,newone,i,jtrue) && 
check(ltrue,newone,i,ktrue)
				&& check(i,newone,jtrue,ktrue) && 
check(ltrue,newone,jtrue,ktrue))
				 temp = temp + 1;
			if (temp == 0) if (check(i,ltrue,jtrue,newone) && 
check(i,ltrue,ktrue,newone)
				&& check(i,newone,jtrue,ktrue) && 
check(ltrue,newone,jtrue,ktrue))
				 temp = temp + 1;
			if (temp == 0)  if (check(i,ltrue,jtrue,newone) && 
check(i,ltrue,ktrue,newone)
				&& check(i,ktrue,jtrue,newone) && 
check(ltrue,ktrue,jtrue,newone))
				 temp = temp + 1;
			if (temp == 0)  if (check(i,ltrue,jtrue,newone) && 
check(i,ltrue,ktrue,newone)
				&& check(i,jtrue,ktrue,newone) && 
check(ltrue,jtrue,ktrue,newone))
				 temp = temp + 1;	
			if (temp == 0)  demerit4 = demerit4 + 1;	
			
			};//	{if newone}

		
		if ((demerit2 < demerit3) && (demerit2 < demerit4)) 
			{
			t4code[i-1][j-2][k-3][l-4][2] = 2;
			totaldemerits = totaldemerits + demerit2;
			}
		else if ((demerit3 < demerit2) && (demerit3 < demerit4)) 
			{
			t4code[i-1][j-2][k-3][l-4][2] = 3;
			totaldemerits = totaldemerits + demerit3;
			}
			
		else if ((demerit4 < demerit2) && (demerit4 < demerit3)) 
			{
			t4code[i-1][j-2][k-3][l-4][2] = 4;
			totaldemerits = totaldemerits + demerit4;
			}
		else
			{
			t4code[i-1][j-2][k-3][l-4][2] = t4code[i-1][j-2][k-3][l-4][1];
			if (t4code[i-1][j-2][k-3][l-4][1]== 2) totaldemerits = totaldemerits + 
demerit2;
			else if (t4code[i-1][j-2][k-3][l-4][1]== 3) totaldemerits = 
totaldemerits + demerit3;
			else totaldemerits = totaldemerits + demerit4;

			};
		

    	
    	};  //{for l = 1 to limit do}

//	{whether or not immediate, now code[1] && code[2] contain successive versions}    
	
 	same = true;
    for (i = 1; i<= limit; i++)
    for (j = i+1; j<= limit; j++)
    for (k = j+1; k<= limit; k++)
    for (l = k+1; l <= limit; l++)		
    	if (same) 	
		if (t4code[i-1][j-2][k-3][l-4][1] != t4code[i-1][j-2][k-3][l-4][2]) same = 
false;

   	
    	
 outcomment <<"Number of demerits in pass " << rep << " was " << totaldemerits << '\n';
  if (same)  outcomment <<"No change in pass " << rep <<" from previous list of 
quartets.\n";
  if (same)  cout <<"No change in pass " << rep <<" from previous list of quartets.\n";
  if (same) actualrepcount = rep;
cout << "Number of demerits in pass " << rep << " was " << totaldemerits << "\n";
 
    for (i = 1; i<= limit; i++)
    for (j = i+1; j<= limit; j++)
    for (k = j+1; k<= limit; k++)
    for (l = k+1; l <= limit; l++)		
		t4code[i-1][j-2][k-3][l-4][1] = t4code[i-1][j-2][k-3][l-4][2];

    	
 } // {for rep = 1 to repcount}   
 
 if (! same) actualrepcount = repcount;
 
 

//{write out modified t4list}

	rep = repcount+1;

    for (i = 1; i<= limit; i++)
    for (j = i+1; j<= limit; j++)
    for (k = j+1; k<= limit; k++)
    for (l = k+1; l <= limit; l++)		
	{
	
	if (t4code[i-1][j-2][k-3][l-4][1] == 2) 
		{
		outt4list.width(2);
		outt4list  << i<< ' ';
		outt4list.width(2);
		outt4list << j << ' ';
		outt4list.width(2);
		outt4list <<  k << ' ';
		outt4list.width(2);
		outt4list << l<< ' ';
		outt4list.precision(7);
		outt4list <<	t4aabb[i-1][j-2][k-3][l-4]<< ' '
		<< t4abab[i-1][j-2][k-3][l-4]<< ' '<< 
			t4abba[i-1][j-2][k-3][l-4]<< ' ' ;
		}
	else if  (t4code[i-1][j-2][k-3][l-4][1] == 3)
		{
		outt4list.width(2);
		outt4list<< i<< ' ';
		outt4list.width(2);
		outt4list << k<< ' ';
		outt4list.width(2);
		outt4list << j<< ' ';
		outt4list.width(2);
		outt4list << l<< ' ';
		outt4list.precision(7);
		outt4list << 
			t4abab[i-1][j-2][k-3][l-4]<< ' '<< 
			t4aabb[i-1][j-2][k-3][l-4]<< ' '<< 
			t4abba[i-1][j-2][k-3][l-4]<< ' ';
		}
	else
		{
		outt4list.width(2);
		outt4list<< i<< ' ';
		outt4list.width(2);
		outt4list << l<< ' ';
		outt4list.width(2);
		outt4list << j<< ' ';
		outt4list.width(2);
		outt4list << k<< ' ';
		outt4list.precision(7);
		outt4list << 
			t4abba[i-1][j-2][k-3][l-4] << ' '<< 
			t4aabb[i-1][j-2][k-3][l-4] << ' '<< 
			t4abab[i-1][j-2][k-3][l-4] << ' ';
		}
	if (t4code[i-1][j-2][k-3][l-4][1] != t4code[i-1][j-2][k-3][l-4][0])
		{
		outt4list <<"  REVERSED " << t4code[i-1][j-2][k-3][l-4][0];
		numberreversed = numberreversed + 1;
		};
	outt4list << '\n';
	
	}; // {for l = 1 to limit do}

outt4list << " 0  0  0  0  0  0  0 \n\n";

outcomment <<"Number of reversed quartets is " << numberreversed << '\n';
cout <<"Number of reversed quartets is " << numberreversed << '\n';
outt4list <<"Number of reversed quartets is " << numberreversed << '\n';
outcomment <<"(If this number is small, the resulting tree is likely to be close to the 
original tree.)";
outt4list <<"Number of error-correcting runs was " << actualrepcount << '\n';
outcomment <<"Number of error-correcting runs was " << actualrepcount << '\n';
outcomment << "Number of demerits in the last pass was " << totaldemerits << '\n';
outt4list << "Number of demerits in the last pass was " << totaldemerits << '\n';
outcomment << "(If this number of demerits is small, the quartets are close to some 
tree.)\n";

  if (same)
  	{
  	outcomment <<"Two successive identical lists were obtained. \n";
  	outt4list <<"Two successive identical lists were obtained. \n";
  	outcomment <<"Further error correction would not be productive.\n";
  	}
  	else 
  	{
  	outcomment <<"No two successive lists were identical.\n";
  	outt4list <<"No two successive lists were identical.\n";
  	outcomment <<"Further error correction might be productive.\n";
  	};
 
   

nearend:  
	finishup();
	return 0;
} //main