// Implementation of the algorithm that is explained in the paper
// "Algorithmic generation of molecular graphs with large Merrifield-Simmons index"
// Stephan Wagner, June 17, 2006

#include <iostream>
#include <fstream>
#include <list>
#include <cln/integer.h>
#include <cln/io.h>
#include <cln/integer_io.h>

using std::ifstream;
using std::ofstream;
using std::cout;
using std::list;
using namespace cln;

// The class for ternary trees.
// A ternary tree is defined by:
// integers sigma, sigma0;
// pointers to the three subtrees, as well as their sizes and positions in the generated lists.
// interval borders for optimality
class TTree {

public:
  cl_I sigma;
  cl_I sigma0;
  TTree * st1;
  TTree * st2;
  TTree * st3;
  int tree_size1;
  int tree_size2;
  int tree_size3;
  int tree_number1;
  int tree_number2;
  int tree_number3;
  cl_I opt_lower_numer;
  cl_I opt_lower_denom;
  cl_I opt_upper_numer;
  cl_I opt_upper_denom;
  void setOpt(cl_I oln, cl_I old, cl_I oun, cl_I oud);
  void showtree(ofstream *file);
  TTree(TTree *t1, TTree *t2, TTree *t3, int ts1, int ts2, int ts3,
int tn1, int tn2, int tn3);
  TTree();
};

// Generates a tree consisting only of the root.
TTree::TTree()
{
  st1 = 0;
  st2 = 0;
  st3 = 0;
  tree_size1 = 0;
  tree_size2 = 0;
  tree_size3 = 0;
  tree_number1 = 1;
  tree_number2 = 1;
  tree_number3 = 1;
  sigma = 2;
  sigma0 = 1;
  opt_lower_numer = 0;
  opt_lower_denom = 1;
  opt_upper_numer = 1;
  opt_upper_denom = 1;
}

// Generates a ternary tree from three given subtrees.
TTree::TTree(TTree *t1, TTree *t2, TTree *t3, int ts1, int ts2, int ts3,
int tn1, int tn2, int tn3)
{
  st1 = t1;
  st2 = t2;
  st3 = t3;
  tree_size1 = ts1;
  tree_size2 = ts2;
  tree_size3 = ts3;
  tree_number1 = tn1;
  tree_number2 = tn2;
  tree_number3 = tn3;
  sigma = t1->sigma * t2->sigma * t3->sigma + 
          t1->sigma0 * t2->sigma0 * t3->sigma0;
  sigma0 = t1->sigma * t2->sigma * t3->sigma;
}

// Writes the definition of a ternary tree to a file, coded as the sizes and list positions of the subtrees.
void TTree::showtree(ofstream *file)
{
  (*file) << "Subtrees: (" << tree_size1 << "," << tree_number1 << ") " <<
          "(" << tree_size2 << "," << tree_number2 << ") " <<
          "(" << tree_size3 << "," << tree_number3 << "); (sigma,sigma_0): (" <<
          sigma << "," << sigma0 << ")\n";
}

// Sets the interval borders for optimality.
void TTree::setOpt(cl_I oln, cl_I old, cl_I oun, cl_I oud)
{
  opt_lower_numer = oln;
  opt_lower_denom = old;
  opt_upper_numer = oun;
  opt_upper_denom = oud;
}

// Checks whether a triple of subtrees can possibly generate an alpha-optimal ternary tree.
bool potentiallyOptimal(TTree *t1, TTree *t2, TTree *t3)
{
 cl_I up1 = t1->opt_lower_denom * t2->opt_lower_numer * t3->opt_lower_numer
   * t1->sigma * t2->sigma0 * t3->sigma0;
 cl_I up2 = t2->opt_lower_denom * t1->opt_lower_numer * t3->opt_lower_numer
   * t2->sigma * t1->sigma0 * t3->sigma0;
 cl_I up3 = t3->opt_lower_denom * t1->opt_lower_numer * t2->opt_lower_numer
   * t3->sigma * t1->sigma0 * t2->sigma0;
 cl_I lo1 = t1->opt_upper_denom * t2->opt_upper_numer * t3->opt_upper_numer
   * t1->sigma * t2->sigma0 * t3->sigma0;
 cl_I lo2 = t2->opt_upper_denom * t1->opt_upper_numer * t3->opt_upper_numer
   * t2->sigma * t1->sigma0 * t3->sigma0;
 cl_I lo3 = t3->opt_upper_denom * t1->opt_upper_numer * t2->opt_upper_numer
   * t3->sigma * t1->sigma0 * t2->sigma0;
 cl_I up;
 cl_I lo;
 if (up1 >= up2)
   {
     if (up2 >= up3) up = up3;
     else up = up2;
   }
 else
   {
     if (up1 >= up3) up = up3;
     else up = up1;
   }
 if (lo1 <= lo2)
   {
     if (lo2 <= lo3) lo = lo3;
     else lo = lo2;
   }
 else
   {
     if (lo1 <= lo3) lo = lo3;
     else lo = lo1;
   }
 if (up * t1->opt_upper_numer * t2->opt_upper_numer * t3->opt_upper_numer >=
     lo * t1->opt_lower_numer * t2->opt_lower_numer * t3->opt_lower_numer)
 return(true);
 else
   return(false);
}

// Compares the Merrifield-Simmons indices of two trees.
bool operator<(const TTree& t1, const TTree& t2)
{
  return(t1.sigma0 < t2.sigma0);
}

// Calculates the upper envelope for a given list of trees.
list<TTree> filter(list<TTree> treelist)
{
  list<TTree> filtered;
  list<TTree>::iterator treelist_iterator = treelist.begin();
  filtered.push_front(treelist.front());
  (filtered.front()).setOpt(0,1,1,1);
  treelist_iterator++;
  while(treelist_iterator != treelist.end())
    {
      TTree t1 = filtered.back();
      TTree t2 = *(treelist_iterator);
      if (t1.opt_lower_numer * (t2.sigma0-t1.sigma0) > t1.opt_lower_denom * (t1.sigma-t2.sigma))
	{
          filtered.pop_back();
          if (filtered.empty())
            {
              filtered.push_front(t2);
              (filtered.front()).setOpt(0,1,1,1);
              treelist_iterator++;
            }
        }
      else if ((t2.sigma0-t1.sigma0) < (t1.sigma-t2.sigma))
        treelist_iterator++;
      else
        {
          if (t1.sigma0 == t2.sigma0)
            t2.setOpt(t1.opt_lower_numer, t1.opt_lower_denom, t1.opt_upper_numer, t1.opt_upper_denom);
          else
            {
              t1.setOpt(t1.opt_lower_numer, t1.opt_lower_denom,t1.sigma-t2.sigma,t2.sigma0-t1.sigma0);
              t2.setOpt(t1.sigma-t2.sigma,t2.sigma0-t1.sigma0,1,1);
            }
          filtered.push_back(t2);
          treelist_iterator++;         
        }
    }
  return(filtered);
}


list<list<TTree> >::iterator iterator_a;
list<list<TTree> >::iterator iterator_b;
list<list<TTree> >::iterator iterator_c;
list<list<TTree> >::iterator iterator_max;
list<TTree>::iterator inner_iterator_a;
list<TTree>::iterator inner_iterator_b;
list<TTree>::iterator inner_iterator_c;
int m,a,b,c,ia,ib,ic,k; 
list<TTree> newlist;

// Implements the main loop of the algorithm.
// treelistlist is the list of lists L(m).
// L(m) is the list of all alpha-optimal ternary trees with 3m+1 vertices.
// The lists are written to a file "lists.txt".
// The lengths of the lists L(m) are written to a file "lengths.txt".
int main()
{
  list<TTree> treelist0;
  treelist0.push_front(TTree::TTree());
  list<list<TTree> > treelistlist;
  treelistlist.push_front(treelist0);
  ofstream lists("lists.txt");
  lists << "The lists L(m):\n\n";
  ofstream listLengths("lengths.txt");
  listLengths << "Lengths of the lists L(m):\n\n" << "1";
  for(m = 1; m < 101; m++)
   {
     newlist.clear();
     iterator_a = treelistlist.begin();
     iterator_max = treelistlist.end();
     iterator_max--;
     for(a = 0; 3*a < m; a++)
       {
	b = a;
        c = m-a-b-1;
        iterator_b = iterator_a;
        iterator_c = iterator_max;
        bool flag = true;
        while(flag)
          {
            ia = 1;
            inner_iterator_a = (*iterator_a).begin();
            while(inner_iterator_a != (*iterator_a).end())
              {
                if (iterator_a == iterator_b)
                  {
                  inner_iterator_b = inner_iterator_a;
                  ib = ia;
                  }
                else
                  {
                  inner_iterator_b = (*iterator_b).begin();
                  ib = 1;
                  }
                while(inner_iterator_b != (*iterator_b).end())
	          {
                    if (iterator_b == iterator_c)
                      {
                      inner_iterator_c = inner_iterator_b;
                      ic = ib;
                      }
                    else
                      {
                      inner_iterator_c = (*iterator_c).begin();
                      ic = 1;
                      }
                    while(inner_iterator_c != (*iterator_c).end())
	              {
                        if (potentiallyOptimal(&(*inner_iterator_a),&(*inner_iterator_b),
					       &(*inner_iterator_c)))
			  newlist.push_front(TTree::TTree(&(*inner_iterator_a),&(*inner_iterator_b),
							  &(*inner_iterator_c),a,b,c,ia,ib,ic));
                        inner_iterator_c++;
                        ic++;
	              }
                    inner_iterator_b++;
                    ib++;
	          }
                inner_iterator_a++;
                ia++;
              }
            if (iterator_b == iterator_c)
              flag = false;
            iterator_b++;
            b++;
            if (iterator_b == iterator_c)
              flag = false;
            iterator_c--;
            c--;
          }
         iterator_a++;
         iterator_max--;
         iterator_max--;
       }
     newlist.sort();
     list<TTree> filteredlist = filter(newlist);     
     treelistlist.push_back(filteredlist); 
     list<TTree> flist = treelistlist.back();
     list<TTree>::iterator filteredlist_iterator = flist.begin();
     lists << "\nList L(" << m << "):\n\n";
     k = 1;
     while(filteredlist_iterator != flist.end())
        {
          lists << "Tree " << k << ". ";
          (*filteredlist_iterator).showtree(&lists);
          filteredlist_iterator++;
          k++;
        }
      listLengths << ", " << filteredlist.size();
    }

  // Determines the optimal trees of size 3m+3.
  // They are written to a file "optlist3.txt".
  // Additionally, the maximal Merrifield-Simmons indices for all sizes are written to a file "maxsigma.txt".
  int curr_opt_size1, curr_opt_size2, curr_opt_num1, curr_opt_num2;
  cl_I curr_opt_sigma;
  iterator_max = treelistlist.begin();
  ofstream optlist3("optlist3.txt");
  ofstream maxsigma("maxsigma.txt");
  optlist3 << "Optimal trees of size 3m+3:\n\n";
  maxsigma << "List of maximal Merrifield-Simmons indices:\n\n";
  for(m = 0; m < 101; m++)
    {
      maxsigma << ((*iterator_max).front()).sigma << "\n";
      maxsigma << (((*iterator_max).back()).sigma + ((*iterator_max).back()).sigma0) << "\n";
     iterator_a = treelistlist.begin();
     iterator_b = iterator_max;
     iterator_max++;
     for(a = 0; 2*a <= m; a++)
       {
	b = m-a;
        inner_iterator_a = (*iterator_a).begin();
        ia = 1;
        while(inner_iterator_a != (*iterator_a).end())
          {
            if (iterator_a == iterator_b)
              {
	        inner_iterator_b = inner_iterator_a;
                ib = ia;
              }
            else
              {
                inner_iterator_b = (*iterator_b).begin();
                ib = 1;
              }
            while(inner_iterator_b != (*iterator_b).end())
	      {
                if(inner_iterator_a->sigma * inner_iterator_b->sigma + 
                   inner_iterator_a->sigma0 * inner_iterator_b->sigma0 > curr_opt_sigma)
		  {
                    curr_opt_size1 = a;
                    curr_opt_size2 = b;
                    curr_opt_num1 = ia;
                    curr_opt_num2 = ib;
                    curr_opt_sigma = inner_iterator_a->sigma * inner_iterator_b->sigma + 
		      inner_iterator_a->sigma0 * inner_iterator_b->sigma0;
                  }
                inner_iterator_b++;
                ib++;
              }
            inner_iterator_a++;
            ia++;
          }
        iterator_a++;
	iterator_b--;
       }
     optlist3 << "m = " << m << ". Subtrees: (" << curr_opt_size1 << "," << curr_opt_num1 << "), (" << curr_opt_size2 << "," << curr_opt_num2 << "); sigma: " << curr_opt_sigma << "\n";
     maxsigma << curr_opt_sigma << "\n";
    }
  listLengths.close();
  lists.close();
  optlist3.close();
  maxsigma.close();
  return(0);
}