/******************************************************************************
Takes sorted lists of kmers and counts matches from lengths 1-k.
Writes the histograms for each pair to DATA/
 *****************************************************************************/

#include <cstdlib>
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <map>
#include <sstream>
#include <math.h>
#include <set>
#include <dirent.h>
#include <sys/types.h>
#include <omp.h>

#include "Coptim1_3Dmat.h"
#include "bitter.h"
#include "amino_acids.h"
#include "util.h"

using namespace std;

const int TAGLENGTH=20;
const bool MCTR_DEBUG=false;
const bool GENEIDRUN=true;
const int CON_MAP_CUTOFF=3; //"cmc" variable.
                            //For a protein with a hit between an HGT-pair to be spared deletion (due to possible HGT), it has to have at least cmc hits among the reference set.  Therefore, a cmc=15 means at least 15 reference orgs have to have hits for the protein to be considered "conserved" (not delete-able). NOTE: Actually this is wrong! This is not what cmc is! It's the number of hits for a reference set organism, but over all the kmers. So this value can be very high.
const int CON_MAP_SET_CUTOFF=3;

const int HGT_MIN_LENGTH=12;
const int MBKL=1.0; //Default=1.0. If tag_length=20, and MBKL=.75, then private variable mbkl will be set to .75*20=15. 15 will then be the minimum length of a match for its score to be passed to the bitter object, which uses these scores to select a right-hand cutoff for the data in the mdist function.
const bool AM=false;
const bool AAM=false;

struct bits_holder
{
  double d[20];
};

struct krow
{
  vector<bool> r;
};

class mctr
{
 private:
  //Variables necessary for all 3 types of run:
  string bfn;                 //base file name: directory holding all the data (e.g. 4ecoli_run1)
  string path;                //full path until base file name
  int tag_length;             //default=20
  int refsetsize;             //0 means no reference genomes.
  string dir;                 //either "" empty string or "scr"
  int filtering_steps;        //=-1 when no filtering used (e.g. using MERGED_TAGS rather than MERGED_TAGS_15_8)
  int filtering_option;       //=-1 when no filtering was used (same as above)
  string merged_directory;    
  string data_directory;
  int sos_cutoff;              //default = 6.5*tag_length
  double max_single_bit_value;  
  int max_bit_sum;             //e.g. tag_length=20, highest bit store=4, 4*20 = 80
  int mbkl;                    //minimum kmer length for bitter to store it, default=tag_length
  vector<string> files;        
  int number_of_genomes;       
  int number_of_genomes_ref;
  int offset;
  int recursion_counter;
  ofstream logstream;

  char array_c[][TAGLENGTH]; //20 is the tag length.
  int *gnums_arr;
  int *gene_ids_arr;
  
  vector<string> array;
  vector<int> gnums;
  vector<int> gene_ids;
  Coptim1_3Dmat* corr_array;
  vector<bitter> bitterv;
  vector<bits_holder> bitsv;
  amino_acids aas;
  util u;
  int run_type; //0=normal run, 1=ref_set_run, 2=final pair run
  bool partial_run;
  char partial_run_lead;
  int get_aa_int[30]; 
  char get_aa_char[20];

  //Variables necessary for HGT correction on subsequent passes through mctr.
  int prob1;                             //ordinal of one of the 2 problematic proteomes in the whole run.
  int prob2;                             //ordinal of one of the 2 problematic proteomes in the whole run.
  vector<krow> rows; //provides transitive information on HGT matches; if the pair shares a long length match, but one member of the pair has matches with the reference set, this applies to the other, so neither protein will be deleted.
  int rows_cols;
  int rows_rows;

  bool pair_run; //true if run is for just a pair with HGT, no reference set (3rd iteration)
  bool msc; //match_score_cutoff; true if the file "match_score_cutoffs.txt" exists, false otherwise.
  bool msc2;
  int genes_found;
  double min_msc;                        //the minimum msc among all the org pairs.
  map<std::pair<int,int>,int> msc_map;   //delete this when subsequent maps are implemented.
  map<std::pair<int,int>,int> con_map;
  vector< set<int> > con_map_set0;
  vector< set<int> > con_map_set1;
  map<std::pair<int,int>,int> hgt_map;
  vector<int> msc_ref_gnums;
  bool distances;
  vector<string> problem_tags;           //vector for the tags from the 2 proteomes that have the HGT problem. merged into the refset later.
  vector<string> ref_tags;
  vector<int> problem_tags_gnums;
  vector<int> problem_tags_gis;
  bool check_composition(string s);
  void get_problem_tags(vector<string> files,string ordinal);
  void get_problem_tags(vector<string> files,string ordinal, set<int> hgt_genes);
  int ppv_index;                         //what line of the problem pair index vector are we on.
  int ref_index;                         //what line of the reference vector are we on.
  vector<std::pair<int, int> > keys;
  vector<string> info_lines0;
  vector<string> info_lines1;
  vector<string> proteins0;
  vector<string> proteins1;
  void setup_it1(ifstream& infostream);
  void setup_it2(ifstream& infostream);

  void read_in_probfiles(ifstream& infostream,vector<string>& files1, vector<string>& files2);

  int calculate_max_bit_sum();
  void open_info_file(ifstream& infostream);
  void create_directories();
  char get_aa_char_fn(int i);
  void set_bitsv();
  void set_get_aa_arrays();

  void set_array_it0(int i);
  void set_array_it1(int i);
  void set_array_it2(int i);

  int set_gh(bool gh[], vector<int>& gnums, int first_row, int last_row);
  void set_score_array(double * score_array,string seq,int count,int gen_ord);

  void find_matches_it0(int col, int first_row, int last_row, string seq);
  void find_matches_it1(int col, int first_row, int last_row, string seq);
  void find_matches_it2(int col, int first_row, int last_row, string seq);
  void increment_bits_it0(int first_row, int last_row,string seq);
  void increment_bits_it1(int first_row, int last_row,string seq);
  void increment_bits_it2(int first_row, int last_row,string seq);
  int get_aa(char aa);
  void write_out_data_it0();
  void write_out_data_it1();
  void write_out_data_it2();
  
  void read_in_files();
  double get_score(int i1,int i2, string seq);
  void populate_bitterv();
  int get_bitterv_index(int i1, int i2);
  int get_sos_cutoff();
  char get_pr_lead();

 public: 
  //constructors
  mctr(string full_path_par,string dir_par, int filtering_steps_par, int filtering_option_par,bool partial_run_par); //Use for the main run when doing entire set.
  mctr(string full_path_par, string dir_par, int filtering_steps_par, int filtering_options_par, int prob1_par, int prob2_par,vector<string>& ref_tags_par); //Use when testing a pair against a reference set.
  mctr(string full_path_par, string dir_par, int filtering_steps_par, int filtering_option_par, int prob1_par, int prob2_par); //Use when running two proteomes, HGT-genes removed, against each other. Final step.

  //run the main algorithm
  void run_it0();    //main run on whole set.
  int  run_it1();    //run pair against each other and the reference set.
  void run_it2();    //run pair against each other with phage/hgt proteins removed.

  void set_problem_pair(int prob1_par,int prob2_par);
};

//Constructor 1. Main run.
//Count all matches between all pairs present in a MERGED_TAGS directory (and also the corresponding MERGED_TAGSscr directory, but has to be called on the data separately using a different object). 
//This could be for a set without filtering (e.g. MERGED_TAGS plain) or with filering (e.g. MERGED_TAGS_15_8). 
//Without filtering: filtering_steps_par=-1, filtering_option_par=-1
mctr::mctr(string full_path_par, string dir_par, int filtering_steps_par, int filtering_option_par,bool partial_run_par)
{
  cout<<"mctr::mctr(string,string,int,int)"<<endl;
  run_type=0;
  partial_run=partial_run_par;
  set_get_aa_arrays();

  //main variables.
  u.extract_paths(full_path_par,path,bfn);
  if(partial_run)
    {
      partial_run_lead=get_pr_lead();
    }
  u.get_rss(path+"/"+bfn+"/"+bfn+"_info.txt",refsetsize);
  u.get_tl(path+"/"+bfn+"/"+bfn+"_info.txt",tag_length);
  mbkl=int(MBKL*(double)tag_length);
  filtering_steps=filtering_steps_par;
  filtering_option=filtering_option_par;
  dir=dir_par;            //dir="" when running on the main set; dir="scr" when running on scrambled set.
  sos_cutoff=get_sos_cutoff(); 
  u.open_ofile_app(logstream,path+"/"+bfn+"/logfile_cm1.txt","logfile_cm1.txt");

  //Set all pair run variables to false.
  ppv_index=ref_index=0;

  string addon="";
  if(filtering_steps!=-1 && filtering_option!=-1)
    {
      addon="_"+u.int_to_string(filtering_steps)+"_"+u.int_to_string(filtering_option);
    }

  merged_directory = "MERGED_TAGS"+dir+addon;
  data_directory = "DATA"+dir+addon;

  ifstream infostream;
  open_info_file(infostream); //sets tag length and number of genomes.
  infostream.close();

  create_directories();
  set_bitsv();  //should have an option: conserved seq bit scores or all bit scores
  max_bit_sum = calculate_max_bit_sum();
  cout<<"max_bit_sum = "<<max_bit_sum<<endl;

  if(max_bit_sum>0 && number_of_genomes>0)
    {
      corr_array = new Coptim1_3Dmat(number_of_genomes,max_bit_sum);
    }
  else
    {
      cout<<"mctr::mctr(string full_path_par, string dir_par, int filtering_steps_par, int filtering_option_par)"<<endl;
      cout<<"max_bit_sum="<<max_bit_sum<<endl;
      cout<<"number_of_genomes="<<number_of_genomes<<endl;
      exit(1);
    }

  populate_bitterv();
  recursion_counter=0;
  read_in_files();
}

//Constructor 2. Pair run against reference.
mctr::mctr(string full_path_par, string dir_par, int filtering_steps_par, int filtering_option_par, int prob1_par, int prob2_par,vector<string>& ref_tags_par)
{
  set_get_aa_arrays();
  run_type=1;
  u.extract_paths(full_path_par,path,bfn);
  u.get_rss(path+"/"+bfn+"/"+bfn+"_info.txt",refsetsize);
  u.get_tl(path+"/"+bfn+"/"+bfn+"_info.txt",tag_length);
  mbkl=int(MBKL*(double)tag_length);

  filtering_steps=filtering_steps_par;
  filtering_option=filtering_option_par;
  dir=dir_par; //do we need this?
  ref_tags=ref_tags_par;
  
  prob1=prob1_par;
  prob2=prob2_par;
  
  //Set all pair run variables to false.
  ppv_index=0;
  ref_index=0;
  sos_cutoff=get_sos_cutoff();

  u.open_ofile_app(logstream,path+"/"+bfn+"/logfile_cm2.txt","logfile_cm2.txt");

  string addon="";
  if(filtering_steps!=-1 && filtering_option!=-1)
    {
      addon="_"+u.int_to_string(filtering_steps)+"_"+u.int_to_string(filtering_option);
    }
  
  merged_directory = "MERGED_TAGS"+dir+addon;
  data_directory = "DATA"+dir+addon;

  ifstream infostream;
  open_info_file(infostream); //sets tag length and number of genomes.
  int bound=number_of_genomes+number_of_genomes_ref; //Need to set the bound variable here, because these two variables will be changed later to reflect that it's a run only against the reference set.
  setup_it1(infostream);        //Reads in the tags (and sorts them) for the 2 problem proteomes. 
  infostream.close();

  set_bitsv();
  max_bit_sum = calculate_max_bit_sum();

  if(max_bit_sum>0)
    {
      int total=refsetsize;
      if(prob1>=refsetsize)
	{
	  total++;
	}
      if(prob2>=refsetsize)
	{
	  total++;
	}

      number_of_genomes_ref=total;
      if(prob1>=bound || prob2>=bound)
	{
	  cout<<"Invalid genome id: prob1="<<prob1<<" prob2="<<prob2<<endl;
	  cout<<"Check HGT_pairs.txt file for genome id out of bounds"<<endl;
	  exit(1);
	}      

      corr_array = new Coptim1_3Dmat(number_of_genomes_ref,max_bit_sum);
    }
  else
    {
      cout<<"Problem in mctr::mctr(string full_path_par, string dir_par, int filtering_steps_par, int filtering_option_par, int prob1_par, int prob2_par,vector<string>& ref_tags)"<<endl;
      cout<<"max_bit_sum="<<max_bit_sum<<endl;
      exit(1);
    }

  populate_bitterv();
  recursion_counter=0;
}

//Constructor 3. Pair run against just itself.
mctr::mctr(string full_path_par, string dir_par, int filtering_steps_par, int filtering_option_par, int prob1_par, int prob2_par)
{
  set_get_aa_arrays();

  if(MCTR_DEBUG)
    {
      cout<<"mctr::mctr(string full_path_par, string dir_par, int filtering_steps_par, int filtering_option_par, int prob1_par, int prob2_par)"<<endl;
      cout<<"Constructor 3 for HGT pair fix"<<endl;
    }

  run_type = 2;

  u.extract_paths(full_path_par,path,bfn);
  u.get_rss(path+"/"+bfn+"/"+bfn+"_info.txt",refsetsize);
  u.get_tl(path+"/"+bfn+"/"+bfn+"_info.txt",tag_length);
  mbkl=int(MBKL*(double)tag_length);
  filtering_steps=filtering_steps_par;
  filtering_option=filtering_option_par;
  dir=dir_par;

  //Set all pair run variables to false.
  prob1=prob1_par;
  prob2=prob2_par;
  ppv_index=0;
  ref_index=0;
  sos_cutoff=get_sos_cutoff();

  u.open_ofile_app(logstream,path+"/"+bfn+"/logfile_cm3.txt","logfile_cm3.txt");

  string addon="";
  if(filtering_steps!=-1 && filtering_option!=-1)
    {
      addon="_"+u.int_to_string(filtering_steps)+"_"+u.int_to_string(filtering_option);
    }
  
  merged_directory = "MERGED_TAGS"+dir+addon;
  data_directory = "DATA"+dir+addon;

  ifstream infostream;
  open_info_file(infostream); //sets tag length and number of genomes.
  setup_it2(infostream);
  infostream.close();

  set_bitsv();
  max_bit_sum = calculate_max_bit_sum();

  number_of_genomes_ref=0;
  number_of_genomes=2;
  corr_array = new Coptim1_3Dmat(2,max_bit_sum);

  populate_bitterv();

  recursion_counter=0;
}

void mctr::run_it0()
{
  if(MCTR_DEBUG)
    {
      cout<<"mctr::run_it0()"<<endl;
    }

  if(!partial_run)
    {
      for(int i=0; i<20; i++)
	{
	  array.clear();
	  gnums.clear();
	  gene_ids.clear();
	  set_array_it0(i);
	  cout<<"First 3 rows of array: "<<endl;
	  //cout<<array.at(0)<<" "<<gnums.at(0)<<" "<<gene_ids.at(0)<<endl;
	  //cout<<array.at(1)<<" "<<gnums.at(1)<<" "<<gene_ids.at(1)<<endl;
	  //cout<<array.at(2)<<" "<<gnums.at(2)<<" "<<gene_ids.at(2)<<endl;
	  cout<<array[0]<<" "<<gnums[0]<<" "<<gene_ids[0]<<endl;
	  cout<<array[1]<<" "<<gnums[1]<<" "<<gene_ids[1]<<endl;
	  cout<<array[2]<<" "<<gnums[2]<<" "<<gene_ids[2]<<endl;
	  cout<<"Last row of array: "<<endl;
	  cout<<array[array.size()-1]<<" "<<gnums[gnums.size()-1]<<" "<<gene_ids[gene_ids.size()-1]<<endl;
	  cout<<"Second to last row of array: "<<endl;
	  cout<<array[array.size()-2]<<endl;
	  cout<<"Last row of array: "<<endl;
	  cout<<array[array.size()-1]<<endl;
	  cout<<"Entering recursive function find_matches_it0(...)..."<<endl;
	  find_matches_it0(0,0,array.size()-1,"");//, aastr, astr);
	  cout<<"Escaped recursive function find_matches(...)!"<<endl;
	}
    }
  else if(partial_run)
    {
      array.clear();
      gnums.clear();
      gene_ids.clear();
      //      set_array_it0(aas.get_aa(partial_run_lead));
      set_array_it0(get_aa_int[(int)partial_run_lead-65]);
      cout<<"First 3 rows of array: "<<endl;
      cout<<array[0]<<" "<<gnums[0]<<" "<<gene_ids[0]<<endl;
      cout<<array[1]<<" "<<gnums[1]<<" "<<gene_ids[1]<<endl;
      cout<<array[2]<<" "<<gnums[2]<<" "<<gene_ids[2]<<endl;
      cout<<"Last row of array: "<<endl;
      cout<<array[array.size()-1]<<" "<<gnums[gnums.size()-1]<<" "<<gene_ids[gene_ids.size()-1]<<endl;
      cout<<"Second to last row of array: "<<endl;
      cout<<array[array.size()-2]<<endl;
      cout<<"Last row of array: "<<endl;
      cout<<array[array.size()-1]<<endl;
      cout<<"Entering recursive function find_matches_it0(...)..."<<endl;
      find_matches_it0(0,0,array.size()-1,"");//, aastr, astr);
      cout<<"Escaped recursive function find_matches(...)!"<<endl;
    }
  
  cout<<"Done with mctr main algorithm."<<endl;

  write_out_data_it0();
}

/*******************
Code for counting matches between problematic (HGT) pairs and a reference set.
 ******************/
int mctr::run_it1()
{ 
  if(MCTR_DEBUG)
    {
      cout<<"bool mctr::run_it1()"<<endl;
    }

  int total=refsetsize;

  if(prob1>=refsetsize)
    {
      total++;
    }
  if(prob2>=refsetsize)
    {
      total++;
    }

  number_of_genomes=number_of_genomes_ref=total;

  if(MCTR_DEBUG)
    {
      cout<<"total="<<total<<endl;
      cout<<"number_of_genomes="<<number_of_genomes<<endl;
    }
  
  for(int i=0; i<20; i++)
    {
      if(MCTR_DEBUG)
	{
	  cout<<i<<"/"<<20<<endl;
	}

      array.clear();
      gnums.clear();
      gene_ids.clear();
      set_array_it1(i);
      
      find_matches_it1(0,0,array.size()-1,"");
    }

  cout<<"Finished match-counting for "<<prob1<<" "<<prob2<<endl;
  cout<<"con_map.size() = "<<con_map.size()<<endl;
  cout<<"con_map_set0.size() = "<<con_map_set0.size()<<endl;
  cout<<"con_map_set1.size() = "<<con_map_set1.size()<<endl;
  cout<<"hgt_map.size() = "<<hgt_map.size()<<endl;
  
#pragma omp critical
  write_out_data_it1();
  
  return genes_found;
}

void mctr::run_it2()
{
  if(MCTR_DEBUG)
    {
      cout<<"void mctr::run_it2()"<<endl;
#pragma omp critical
      logstream<<"//In void mctr::run_it2() for prob1="<<prob1<<" prob2="<<prob2<<endl;
    }
  
  for(int i=0; i<20; i++)
    {
      array.clear();
      gnums.clear();
      gene_ids.clear();

      set_array_it2(i);

      if(MCTR_DEBUG)
	{
	  cout<<prob1<<" "<<prob2<<" Last row of array: "<<endl;
	  cout<<array[array.size()-1]<<" "<<gnums[gnums.size()-1]<<" "<<gene_ids[gene_ids.size()-1]<<endl;
	  cout<<"Second to last row of array: "<<endl;
	  cout<<array[array.size()-2]<<" "<<gnums[gnums.size()-2]<<" "<<gene_ids[gene_ids.size()-2]<<endl;
	  cout<<prob1<<" "<<prob2<<" "<<msc<<endl;
	}

      find_matches_it2(0,0,array.size()-1,"");

      if(MCTR_DEBUG)
	{
	  cout<<prob1<<" "<<prob2<<" Done in find_matches"<<endl;
	}
     }  

  cout<<"Finished match-counting (it2) for "<<prob1<<" "<<prob2<<endl;

  write_out_data_it2();
}



/*****************************************************
set_bitsv()
*****************************************************/
void mctr::set_bitsv()
{
  logstream<<"void mctr::set_bitsv()"<<endl;
  logstream<<"Setting the bitsv vector..."<<endl;

  if(run_type==0) //the normal case.
    {
      double max=-1;	  
	  
      for(int numg=0; numg<number_of_genomes; numg++)
	{
	  ifstream histstream;
	  u.open_ifile(histstream,path+"/"+bfn+"/HIST/"+bfn+"_"+u.int_to_string(numg)+"_hist.txt","histstream");
	  
	  bits_holder bh_temp;
	  int gnum;
	  histstream>>gnum;
	  if(gnum!=numg)
	    {
	      cout<<"Mismatch in set_bitsv case 0."<<endl;
	      exit(1);
	    }

	  int total_residues;
	  histstream>>total_residues;

	  for(int j=0; j<20; j++)
	    {
	      int count;
	      double frequency;
	      histstream>>count>>frequency;
	      bh_temp.d[j]=(double)count/(double)total_residues;
	      bh_temp.d[j]=( (log(bh_temp.d[j])) * (-1.0)) / 2.0; //NOTE: using log, which in C++ uses base e.
	      //NOTE: Divide by 2 here to speed up calculating scores in get_score(...)
	      if(max==-1 || max<bh_temp.d[j])
		{
		  max = bh_temp.d[j];
		}
	    }
	    bitsv.push_back(bh_temp);
	    histstream.close();
	}
      
      max_single_bit_value=max;
    } //end of normal case run_type==0

  else if(run_type==1) //Running prospective HGT pair against a reference set.
    {
      ofstream outstream;
      u.open_ofile(outstream,path+"/"+bfn+"/"+bfn+"_msc_ref_keys.txt","outstream");

      //Original ordinal, then new ordinal
      outstream<<prob1<<" "<<0<<endl;
      outstream<<prob2<<" "<<1<<endl;
      keys.push_back(std::make_pair(prob1,0) );
      keys.push_back(std::make_pair(prob2,1) );

      for(int i=0; i<refsetsize; i++)
	{
	  if(i!=prob1 && i!=prob2)
	    {
	      int p1,p2;
	      if(i>prob1 && i>prob2)
		{
		  outstream<<i<<" "<<i<<endl;
		  p1=i;
		  p2=i;
		}
	      else if(i>prob1)
		{
		  outstream<<i<<" "<<i+1<<endl;
		  p1=i;
		  p2=i+1;
		}
	      else
		{
		  outstream<<i<<" "<<i+2<<endl;
		  p1=i;
		  p2=i+2;
		}
	      keys.push_back(std::make_pair(p1,p2) );
	    }
	}

      outstream.close();
      
      //get the bits for the 2 problem proteomes.
      double max=-1;
      for(int i=0; i<number_of_genomes; i++) //number_of_genomes, NOT number_of_genomes_ref, because the code has to go through the entire set to find a) all the reference genomes and b) the two problem genomes.
	{
	  ifstream probstream;
	  string probfile=path+"/"+bfn+"/HIST/"+bfn+"_"+u.int_to_string(i)+"_hist.txt";
	  probstream.open(probfile.c_str());
	  if(probstream.fail())
	    {
	      cout<<"Failed to open probstream to probfile="<<probfile<<endl;
	      cout<<"run_type=="<<run_type<<endl;
	      exit(1);
	    }
	  
	  int gnum;
	  probstream>>gnum;
	  int total_residues;
	  probstream>>total_residues;

	  bits_holder bh_temp;
	  for(int j=0; j<20; j++)
	    {
	      int count;
	      double frequency;
	      probstream>>count>>frequency;
	      bh_temp.d[j]=(double)count/(double)total_residues;
	      bh_temp.d[j]=( (log(bh_temp.d[j])) * (-1.0)) / 2.0; //divide by 2 here to speed up calculating scores in get_score(...)
	      if(max==-1 || max<bh_temp.d[j])
		{
		  max = bh_temp.d[j];
		}
	    }

	  if(gnum==prob1)
	    {
	      bitsv.insert(bitsv.begin(),bh_temp);
	    }
	  else if(gnum==prob2)
	    {
	      vector<bits_holder>::iterator it;
	      it=bitsv.begin();
	      it++;
	      bitsv.insert(it,bh_temp);
	    }
	  else if(gnum<refsetsize)
	    {
	      bitsv.push_back(bh_temp);
	    }
	  probstream.close();      
	}

      max_single_bit_value=max;
    } //run_type==1

  else if(run_type==2) //Running a pair against itself, with the suspicious genes already flagged by the previous run against a reference set.
    {
      //get the bits for the 2 problem proteomes.
      double max=-1;
      for(int i=0; i<number_of_genomes; i++) //number_of_genomes, NOT number_of_genomes_ref, because the code has to go through the entire set to find a) all the reference genomes and b) the two problem genomes.
	{
	  ifstream probstream;
	  string probfile=path+"/"+bfn+"/HIST/"+bfn+"_"+u.int_to_string(i)+"_hist.txt";
	  probstream.open(probfile.c_str());
	  if(probstream.fail())
	    {
	      cout<<"Failed to open probstream to probfile="<<probfile<<endl;
	      cout<<"run_type=="<<run_type<<endl;
	      exit(1);
	    }

	  bits_holder bh_temp;
	  int gnum;
	  probstream>>gnum;
	  int total_residues;
	  probstream>>total_residues;

	  for(int j=0; j<20; j++)
	    {
	      int count;
	      double frequency;
	      probstream>>count>>frequency;
	      bh_temp.d[j]=(double)count/(double)total_residues;
	      bh_temp.d[j]=( (log(bh_temp.d[j])) * (-1.0)) / 2.0; //divide by 2 here to speed up calculating scores in get_score(...)
	      if(max==-1 || max<bh_temp.d[j])
		{
		  max = bh_temp.d[j];
		}
	    }
	  
	  if(gnum==prob1)
	    {
	      bitsv.insert(bitsv.begin(),bh_temp);
	    }
	  else if(gnum==prob2)
	    {
	      vector<bits_holder>::iterator it;
	      it=bitsv.begin();
	      it++;
	      bitsv.insert(it,bh_temp);
	    }
	  probstream.close();
	}
      
      max_single_bit_value=max;
      
      if(bitsv.size()!=2)
	{
	  cout<<"Problem in void mctr::set_bitsv()"<<endl;
	  cout<<"bitsv.size()="<<bitsv.size()<<endl;
	  cout<<"prob1 = "<<prob1<<endl;
	  cout<<"prob2 = "<<prob2<<endl;
	  exit(1);
	}
    } //run_type==2

} //end of set_bitsv()

/***********************************
set_array()
************************************/
void mctr::set_array_it0(int i)
{
  cout<<"In mctr_optim1::set_array(int i). int i = "<<i<<endl;
  logstream<<"In mctr_optim1::set_array(int i). int i = "<<i<<endl;
  //char current_amino_acid = aas.get_aa(i);
  stringstream ss;
  ss<<i;
  string s = ss.str();
  ifstream instream;
  string infile = path + "/"+bfn + "/"+merged_directory+"/TAGS"+s+".txt";
  cout<<"infile = "<<infile<<endl;
  logstream<<"infile = "<<infile<<endl;
  instream.open( infile.c_str());
  if(instream.fail())
    {
      cout<<"In mctr::set_array(int i)"<<endl;
      cout<<"Failed to open instream to infile = "<<infile<<endl;
      exit(1);
    }
  
  int total_size;
  instream>>total_size;
  cout<<"\tTotal tags in current file = "<<total_size<<endl;
  logstream<<"\tTotal tags in current file = "<<total_size<<endl;
  
  for(int k=0; k<total_size; k++)
    {
      string sequence;
      int gnumber;
      int gene_id;
      instream>>sequence>>gnumber;
      if(sequence.length()<tag_length)
	{
	  cout<<"Hit the empty string. Bye."<<endl;
	  exit(1);
	}
      instream>>gene_id;
      gene_ids.push_back(gene_id);	   
      array.push_back(sequence);
      gnums.push_back(gnumber);
    }
  
  instream.close();
  cout<<"Closed "<<infile<<endl;
  cout<<"Size of array = "<<array.size()<<endl;
  cout<<"Size of gnums = "<<gnums.size()<<endl;      
  cout<<"Size of gene_ids = "<<gene_ids.size()<<endl;
  
  if(array.size()!=gnums.size())
    {
      cout<<"array.size()!=gnums.size()"<<endl;
      exit(1);
    }
  
  else if(gene_ids.size()!=gnums.size())
    {
      cout<<"gene_ids.size()="<<gene_ids.size()<<" != gnums.size()="<<gnums.size()<<endl;
      exit(1);
    }
} // *****end of mctr::set_array_it0(int i)*****

/***********************************
set_array_it1()
************************************/
void mctr::set_array_it1(int i)
{
  if(MCTR_DEBUG)
    {
      cout<<"void mctr::set_array_it1(int i)"<<endl;
      cout<<"i="<<i<<endl;
      cout<<"ref_index="<<ref_index<<endl;
      cout<<"ppv_index="<<ppv_index<<endl;
      cout<<"ref_tags.size()="<<ref_tags.size()<<endl;
    }
  
  logstream<<"void mctr::set_array_it1(int i)"<<endl;
  logstream<<"i="<<i<<endl;

  //  char caa = aas.get_aa(i);
  char caa = get_aa_char[i];

  string rs=ref_tags[ref_index]; //set to 0 in constructor.
  string ps=problem_tags[ppv_index];
  
  if(MCTR_DEBUG)
    {
      cout<<"caa="<<caa<<endl;
      cout<<"rs="<<rs<<endl;
      cout<<"ps="<<ps<<endl;
    }
  
  while(rs[0]==caa && ps[0]==caa && ref_index<ref_tags.size() && ppv_index<problem_tags.size() )
    {
      if(rs<ps)
	{
	  string sequence;
	  int gnumber;
	  int gene_id;
	  istringstream iss(rs);
	  iss>>sequence>>gnumber>>gene_id;
	  //leave out tags from problem sequences since these problem proteomes have already had their tags generated in problem_tags.
	  int new_number = gnumber-2;
	  if( new_number>prob2)
	    {
	      //new_number=gnumber-2;
	    }
	  else if(new_number==prob2)
	    {
	      new_number=-1;
	    }
	  else if(new_number>prob1)
	    {
	      new_number=gnumber-1;
	    }
	  else if(new_number==prob1)
	    {
	      new_number=-1;
	    }
	  else if(new_number<prob1)
	    {
	      new_number=gnumber;
	    }
	  
	  if(new_number!=-1)
	    {
	      array.push_back(sequence);
	      gnums.push_back(new_number);
	      gene_ids.push_back(gene_id);
	    }
	  
	  ref_index++;
	  if(ref_index<ref_tags.size())
	    {
	      //	      rs=ref_tags.at(ref_index);
	      rs=ref_tags[ref_index];
	    }
	}
      else
	{
	  string sequence;
	  int gnumber;
	  int gene_id;
	  istringstream iss(ps);
	  iss>>sequence>>gnumber>>gene_id;
	  array.push_back(sequence);
	  gnums.push_back(gnumber);
	  gene_ids.push_back(gene_id);
	  ppv_index++;
	  if(ppv_index<problem_tags.size())
	    {
	      //	      ps=problem_tags.at(ppv_index);
	      ps=problem_tags[ppv_index];
	    }
	}
    }
  while(rs[0]==caa && ref_index<ref_tags.size() )
    {
      string sequence;
      int gnumber;
      int gene_id;
      istringstream iss(rs);
      iss>>sequence>>gnumber>>gene_id;
      int new_number = gnumber-2;
      if( new_number>prob2)
	{
	  //new_number=gnumber-2;
	}
      else if(new_number==prob2)
	{
	  new_number=-1;
	}
      else if(new_number>prob1)
	{
	  new_number=gnumber-1;
	}
      else if(new_number==prob1)
	{
	  new_number=-1;
	}
      else if(new_number<prob1)
	{
	  new_number=gnumber;
	}
      
      if(new_number!=-1)
	{
	  array.push_back(sequence);
	  gnums.push_back(new_number);
	  gene_ids.push_back(gene_id);
	}
      
      ref_index++;
      if(ref_index<ref_tags.size())
	{
	  //	  rs=ref_tags.at(ref_index);
	  rs=ref_tags[ref_index];
	}
    }
  while(ps[0]==caa && ppv_index<problem_tags.size())
    {
      string sequence;
      int gnumber;
      int gene_id;
      istringstream iss(ps);
      iss>>sequence>>gnumber>>gene_id;
      array.push_back(sequence);
      gnums.push_back(gnumber);
      gene_ids.push_back(gene_id);
      ppv_index++;
      if(ppv_index<problem_tags.size())
	{
	  //	  ps=problem_tags.at(ppv_index);
	  ps=problem_tags[ppv_index];
	}
    }
  
  cout<<prob1<<" "<<prob2<<": array set ("<<caa<<")."<<endl;

  if(array.size()!=gnums.size())
    {
      cout<<"array.size()!=gnums.size()"<<endl;
      exit(1);
    }
  
  else if(gene_ids.size()!=gnums.size())
    {
      cout<<"gene_ids.size()="<<gene_ids.size()<<" != gnums.size()="<<gnums.size()<<endl;
      exit(1);
    }
} // *****end of mctr::set_array_it1(int i)*****

/***********************************
set_array_2()
************************************/
void mctr::set_array_it2(int i)
{
  logstream<<"In mctr_optim1::set_array(int i). int i = "<<i<<endl;

  //  char caa = aas.get_aa(i);
  char caa=get_aa_char[i];

  //  string ps=problem_tags.at(ppv_index);
  string ps=problem_tags[ppv_index];
  
  while(ps[0]==caa && ppv_index<problem_tags.size())
    {
      string sequence;
      int gnumber;
      int gene_id;
      istringstream iss(ps);
      iss>>sequence>>gnumber>>gene_id;
      array.push_back(sequence);
      gnums.push_back(gnumber);
      gene_ids.push_back(gene_id);
      ppv_index++;
      if(ppv_index<problem_tags.size())
	{
	  //	  ps=problem_tags.at(ppv_index);
	  ps=problem_tags[ppv_index];
	}
    }

  cout<<prob1<<" "<<prob2<<": array set(it2)"<<endl;

} // *****end of mctr::set_array_it2(int i)*****


/************************************************
find_matches_it0(...)
Count matches between all pairs. 
No correction for HGT at this stage: the original match-counting.
**************************************************/
void mctr::find_matches_it0(int col, int first_row, int last_row, string seq)
{
  if(MCTR_DEBUG)
    {
      cout<<"find_matches_it0(...): col = "<<col<<" first_row = "<<first_row<<" last_row = "<<last_row<<" seq = "<<seq<<endl;
      //cout<<"array.at(first_row) = " <<array.at(first_row)<<endl;
      cout<<"array.at(first_row) = " <<array[first_row]<<endl;
      //      cout<<"array.at(last_row) = "<<array.at(last_row)<<endl;
      cout<<"array.at(last_row) = "<<array[last_row]<<endl;
    }

  recursion_counter++;
  if(recursion_counter%100000==0)
    {
      cout<<"first_row = "<<first_row<<endl;
      //      cout<<array.at(first_row)<<" "<<gnums.at(first_row)<<" "<<gene_ids.at(first_row)<<endl;
      cout<<array[first_row]<<" "<<gnums[first_row]<<" "<<gene_ids[first_row]<<endl;
    }

  //IF WE'VE ONLY GOT ONE SEQUENCE
  if(first_row==last_row)
    {
      //do nothing
    }

  //MULTIPLE SEQUENCES, LAST COLUMN
  else if(col==tag_length)
    {
      //do nothing, past the end of the sequence (e.g. previous recursion hit 20)
    }

  //IF WE'VE GOT MULTIPLE SEQUENCES (I.E. MORE MATCHES), NOT MAX LENGTH
  else
    {
      //SHORTCUT: IF THE ENTIRE BLOCK HAS THE SAME LEADING CHARACTER
      //      if(array.at(last_row)[col]==array.at(first_row)[col])
      if(array[last_row][col]==array[first_row][col])
	{
	  //	  int aa_index=get_aa(array.at(first_row)[col]);
	  //	  int aa_index=get_aa_int[(int)(array.at(first_row)[col])-65];
	  int aa_index=get_aa_int[(int)(array[first_row][col])-65];
	  
	  //IF THE SINGLE-CHARACTER COL IS NOT A ^ COLUMN
	  if(aa_index>=0 && aa_index<20)
	    {
	      find_matches_it0(col+1,first_row,last_row,seq+array[first_row][col]);
	      increment_bits_it0(first_row,last_row,seq+array[first_row][col]);
	    }
	}

      //IF THERE ARE DIFFERENT LEADING CHARACTERS, NEED TO PARTITION ARRAY FURTHER
      else
	{
	  int index_start=first_row;
	  int index_end=first_row;
	  
	  while(index_start<last_row)
	    {
	      //	      char leader=array.at(index_start)[col];
	      char leader=array[index_start][col];
	      
	      while(index_end<=last_row && array[index_end][col]==leader)
		{
		  index_end++;
		}
	      
	      //	      int aa_index = get_aa(leader);
	      int aa_index = get_aa_int[(int)leader-65];
	      
	      //IF THEY'RE THE SAME VALUE (SHOULD NEVER HAPPEN)
	      if(index_end==index_start)
		{
		}
	      //IF THEY'RE JUST 1 OFF: NO MATCHES AT THIS DEPTH
	      else if(index_end==index_start+1)
		{
		}
	      //IF THERE ARE MATCHES
	      else if(aa_index<20 && aa_index>=0)
		{
		  find_matches_it0(col+1, index_start, index_end-1,seq+leader);
		  increment_bits_it0(index_start,index_end-1,seq+leader);
		}
	      index_start=index_end;	      
	    } //end of while(index_start<last_row)
	}
    }
  if(MCTR_DEBUG)
    {
      cout<<"End of find_matches_it0(...)"<<endl;
    }
}

/************************************************
find_matches_it1(...)
Run between problematic (HGT) pair and reference set.
Match-counting not really used.
Just quickly assessing and marking long kmer matches in the set.
**************************************************/
void mctr::find_matches_it1(int col, int first_row, int last_row, string seq)
{
  if(MCTR_DEBUG)
    {
      cout<<"col = "<<col<<" first_row = "<<first_row<<" last_row = "<<last_row<<" seq = "<<seq<<endl;
      //      cout<<"array.at(first_row) = " <<array.at(first_row)<<endl;
      //      cout<<"array.at(last_row) = "<<array.at(last_row)<<endl;
      cout<<"array.at(first_row) = " <<array[first_row]<<endl;
      cout<<"array.at(last_row) = "<<array[last_row]<<endl;
    }

  recursion_counter++;

  //IF WE'VE ONLY GOT ONE SEQUENCE
  if(first_row==last_row)
    {
      //do nothing
    }

  //MULTIPLE SEQUENCES, LAST COLUMN
  else if(col==tag_length)
    {
      //do nothing, past the end of the sequence (e.g. previous recursion hit 20)
    }

  //IF WE'VE GOT MULTIPLE SEQUENCES (I.E. MORE MATCHES), NOT MAX LENGTH
  else
    {
      //SHORTCUT: IF THE ENTIRE BLOCK HAS THE SAME LEADING CHARACTER
      //      if(array.at(last_row)[col]==array.at(first_row)[col])
      if(array[last_row][col]==array[first_row][col])
	{
	  //	  int aa_index=get_aa(array.at(first_row)[col]);
	  //	  int aa_index=get_aa_int[(int)(array.at(first_row)[col])-65];
	  int aa_index=get_aa_int[(int)(array[first_row][col])-65];
	  
	  //IF THE SINGLE-CHARACTER COL IS NOT A ^ COLUMN
	  if(aa_index>=0 && aa_index<20)
	    {
	      //	      find_matches_it1(col+1,first_row,last_row,seq+array.at(first_row)[col]);
	      //	      increment_bits_it1(first_row,last_row,seq+array.at(first_row)[col]);
	      find_matches_it1(col+1,first_row,last_row,seq+array[first_row][col]);
	      increment_bits_it1(first_row,last_row,seq+array[first_row][col]);
	    }
	}//end of shortcut
      
      //IF THERE ARE DIFFERENT LEADING CHARACTERS, NEED TO PARTITION ARRAY FURTHER
      else
	{
	  int index_start=first_row;
	  int index_end=first_row;
	  
	  while(index_start<last_row)
	    {
	      //	      char leader=array.at(index_start)[col];
	      char leader=array[index_start][col];
	      
	      while(index_end<=last_row && array[index_end][col]==leader)
		{
		  index_end++;
		}
	      	      
	      //	      int aa_index = get_aa(leader);
	      int aa_index = get_aa_int[(int)leader-65];

	      //IF THEY'RE THE SAME VALUE (SHOULD NEVER HAPPEN)
	      if(index_end==index_start)
		{
		}
	      //IF THEY'RE JUST 1 OFF: NO MATCHES AT THIS DEPTH
	      else if(index_end==index_start+1)
		{
		}
	      //IF THERE ARE MATCHES
	      else if(aa_index<20 && aa_index>=0)
		{
		  find_matches_it1(col+1, index_start, index_end-1,seq+leader);
		  increment_bits_it1(index_start,index_end-1,seq+leader);
		}
	      index_start=index_end;	      
	    } //end of while(index_start<last_row)
	}
    }
} //end of find_matches_it1(...)

/************************************************
find_matches_it2(...)
**************************************************/
void mctr::find_matches_it2(int col, int first_row, int last_row, string seq)//, ofstream& alalstream, ofstream& alstream)
{
  if(MCTR_DEBUG)
    {
      cout<<"col = "<<col<<" first_row = "<<first_row<<" last_row = "<<last_row<<" seq = "<<seq<<endl;
      //      cout<<"array.at(first_row) = " <<array.at(first_row)<<endl;
      //      cout<<"array.at(last_row) = "<<array.at(last_row)<<endl;
      cout<<"array.at(first_row) = " <<array[first_row]<<endl;
      cout<<"array.at(last_row) = "<<array[last_row]<<endl;
    }

  recursion_counter++;
  if(recursion_counter%100000==0)
    {
      cout<<"first_row = "<<first_row<<" for prob1="<<prob1<<" prob2="<<prob2<<endl;
      //      cout<<array.at(first_row)<<" "<<gnums.at(first_row)<<" "<<gene_ids.at(first_row)<<endl;
      cout<<array[first_row]<<" "<<gnums[first_row]<<" "<<gene_ids[first_row]<<endl;
    }

  //IF WE'VE ONLY GOT ONE SEQUENCE
  if(first_row==last_row)
    {
      //do nothing
    }

  //MULTIPLE SEQUENCES, LAST COLUMN
  else if(col==tag_length)
    {
      //do nothing, past the end of the sequence (e.g. previous recursion hit 20)
    }

  //IF WE'VE GOT MULTIPLE SEQUENCES (I.E. MORE MATCHES), NOT MAX LENGTH
  else
    {
      //SHORTCUT: IF THE ENTIRE BLOCK HAS THE SAME LEADING CHARACTER
      //      if(array.at(last_row)[col]==array.at(first_row)[col])
      if(array[last_row][col]==array[first_row][col])
	{
	  //	  int aa_index=get_aa(array.at(first_row)[col]);
	  //	  int aa_index=get_aa_int[(int)array.at(first_row)[col]-65];
	  int aa_index=get_aa_int[(int)array[first_row][col]-65];
	  
	  //IF THE SINGLE-CHARACTER COL IS NOT A ^ COLUMN
	  if(aa_index>=0 && aa_index<20)
	    {
	      //	      find_matches_it2(col+1,first_row,last_row,seq+array.at(first_row)[col]);
	      //	      increment_bits_it2(first_row,last_row,seq+array.at(first_row)[col]);
	      find_matches_it2(col+1,first_row,last_row,seq+array[first_row][col]);
	      increment_bits_it2(first_row,last_row,seq+array[first_row][col]);
	    }
	}//end of shortcut

      //IF THERE ARE DIFFERENT LEADING CHARACTERS, NEED TO PARTITION ARRAY FURTHER
      else
	{
	  int index_start=first_row;
	  int index_end=first_row;

	  while(index_start<last_row)
	    {
	      char leader=array[index_start][col];
	      
	      while(index_end<=last_row && array[index_end][col]==leader)
		{
		  index_end++;
		}
	      	      
	      //	      int aa_index = get_aa(leader);
	      int aa_index = get_aa_int[(int)leader-65];

	      //IF THEY'RE THE SAME VALUE (SHOULD NEVER HAPPEN)
	      if(index_end==index_start)
		{
		}
	      //IF THEY'RE JUST 1 OFF: NO MATCHES AT THIS DEPTH
	      else if(index_end==index_start+1)
		{
		}
	      //IF THERE ARE MATCHES
	      else if(aa_index<20 && aa_index>=0)
		{
		  find_matches_it2(col+1, index_start, index_end-1,seq+leader);
		  increment_bits_it2(index_start,index_end-1,seq+leader);
		}
	      index_start=index_end;	      
	    } //end of while(index_start<last_row)
	}
    }
}

/*******************************************
increment_bits_it0(...)
-used for run_it0; the main run, counting all matches between all inputs.
-separate from hgt/phage fix.
*********************************************/
void mctr::increment_bits_it0(int first_row, int last_row,string seq)//,ofstream& alalstream, ofstream& alstream)
{
  if(MCTR_DEBUG)
    {
      cout<<"increment_bits_it0(...)"<<endl;
      cout<<"first_row="<<first_row<<" last_row="<<last_row<<" seq="<<seq<<endl;
    }

  bool * gh; //replacing gnum_holder.
  gh = new bool[number_of_genomes];
  //  cout<<"number of genomes="<<number_of_genomes<<endl;
  
  for(int i=0; i<number_of_genomes; i++)
    {
      gh[i]=false;
    }

  int num_uniq_genomes = set_gh(gh,gnums,first_row,last_row);

  //  cout<<"Top of increment_bits_it0..."<<endl;
  //  cout<<"size of gh = number_of_genomes = "<<number_of_genomes<<endl;
  //  cout<<"num_uniq_genomes = "<<num_uniq_genomes<<endl;
  //  cout<<"Contents of gh:"<<endl;
  //  for(int i=0; i<number_of_genomes; i++)
  //   {
  ///      cout<<i<<" "<<gh[i]<<endl;
  //    }
  
  //  set<int> gnum_holder; //using a set to deal with repeats within the same organism (sets are non-redundant)
  //bool store_gene_ids=false; //Used for storing gene numbers in bitterv because their score >= match_score_cutoff.

  //  bool full=false;
  //set<int> gnum_diag_hit;
  //pair<set<int>::iterator,bool> ret;

  //Store all ordinal IDs in a set.
  /*  for(int i=first_row; i<=last_row && !full; i++)
    {
      //      int gnum=gnums.at(i);
      int gnum=gnums[i];
      ret = gnum_holder.insert(gnum);
      if(gnum_holder.size()>=number_of_genomes && gnum_diag_hit.size()>=number_of_genomes)
	{
	  full=true;
	}
      if(ret.second==false) //This means gnum was already in the set. Since the sequence is present twice in the same organism, that means the organism has a match within itself; therefore add the number to the diagonal set.
	{
	  gnum_diag_hit.insert(gnum);
	}
    }
  */

  if(MCTR_DEBUG)
    {
      cout<<"In increment_bits_it0(...)"<<endl;
      cout<<"first_row = "<<first_row<<" last_row = "<<last_row<<" seq = "<<seq<<endl;
      //     cout<<"ret = ";
      //set<int>::iterator it;
      /*      for(it=gnum_holder.begin(); it!=gnum_holder.end(); it++)
	{
	  cout<<*it<<" ";
	}
      */
      //cout<<endl;
    }

  //Case 1: No matches between organisms or within organisms.
  //  if(gnum_holder.size()<=1 && gnum_diag_hit.size()==0)
  if(num_uniq_genomes<=1)
    {
      //      cout<<"Doing nothing"<<endl;
      //do nothing
    }
  //Case 2: A match within an organism and not between organisms
  /*  else if(gnum_holder.size()<=1 && gnum_diag_hit.size()>=1)
    {
      set<int>::iterator it3;
      for(it3=gnum_diag_hit.begin(); it3!=gnum_diag_hit.end(); it3++)
	{
	  int it3_int = *it3;
	  double score = get_score(it3_int,it3_int,seq);
	  corr_array->increment_element(it3_int,it3_int,(int)(score+.5));
	}
	}*/ 
  
  //Case 3: A match between ALL organisms
  //  else if(gnum_holder.size()>=number_of_genomes) //shortcut
  else if(num_uniq_genomes==number_of_genomes)
    {
      //      cout<<"full sized"<<endl;
      double * score_array;
      score_array = new double[seq.length()*num_uniq_genomes];
      
      //set a vector of scores
      int count=0;	  
      for(int i=0; i<number_of_genomes; i++)
	{
	  if(gh[i])
	    {
	      set_score_array(score_array,seq,count,i);
	      count++;
	    }
	}
      
      for(int i=0; i<number_of_genomes-1; i++)
	{
	  for(int j=i+1; j<number_of_genomes; j++)
	    {
	      double score=0;
	      
	      for(int k=0; k<seq.length(); k++)
		{
		  score+=score_array[(i*seq.length())+k]+score_array[(j*seq.length())+k];
		}
	      
	      corr_array->increment_element(i,j,(int)(score+.5));
	      if(seq.length()>=mbkl)
		{
		  bitterv[get_bitterv_index(i,j)].add_new_score(score);
		}
	    }
	}
      
      /*      for(int i=0; i<number_of_genomes-1; i++)
	{	  
	  for(int j=i+1; j<number_of_genomes; j++)
	    {
	      double score = get_score(i,j,seq);
	      corr_array->increment_element(i,j,(int)(score+.5));
	      if(seq.length()>=mbkl)
		{
//		  bitterv.at(get_bitterv_index(i,j)).add_new_score(score);
		  bitterv[get_bitterv_index(i,j)].add_new_score(score);
		}
	    }
	}
      */
	/*      if(gnum_diag_hit.size()>0)
	{
	  set<int>::iterator it1;
	  for(it1=gnum_diag_hit.begin(); it1!=gnum_diag_hit.end(); it1++)
	    {
	      int it1_int = *it1;
	      double score = get_score(it1_int,it1_int,seq);
	      corr_array->increment_element(it1_int,it1_int,(int)(score+.5));
	    }
	    }*/
      delete[] score_array;
    }

  //Case 4: A match between SOME non-identical organisms
  else
    {
      /*      cout<<"Case 4: Some but not all have matches."<<endl;
      cout<<"num_uniq_genomes = "<<num_uniq_genomes<<endl;
      cout<<"number_of_genomes = "<<number_of_genomes<<endl;
      cout<<"seq="<<seq<<endl;
      */
      
      //new code
      double * score_array;
      score_array = new double[seq.length()*num_uniq_genomes];
      //      cout<<"Size of score_array = "<<seq.length()*num_uniq_genomes<<endl;
      
      //set a vector of scores
      //      for(int i=0; i<number_of_genomes; i++)
      //    cout<<"In Case 4: going through set_score_array loop..."<<endl;
      int count=0;
      for(int i=0; i<number_of_genomes; i++)
	{
	  //	  cout<<"i="<<i<<endl;
	  if(gh[i])
	    {
	      set_score_array(score_array,seq,count,i);
	      count++;
	    }
	}
      //      cout<<"Done in loop"<<endl;
      
      int gh_index1=0;
      //      cout<<"new loop..."<<endl;

      for(int i=0; i<number_of_genomes-1; i++)
	{
	  //	  cout<<"i="<<i<<" out of "<<number_of_genomes-1<<endl;
	  if(gh[i])
	    {
	      int gh_index2=gh_index1+1;
	      for(int j=i+1; j<number_of_genomes; j++)
		{
		  //		  cout<<"j="<<j<<" out of "<<number_of_genomes<<endl;
		  if(gh[j])
		    {
		      double score=0;
		      for(int k=0; k<seq.length(); k++)
			{
			  score+=score_array[(gh_index1*seq.length())+k]+score_array[(gh_index2*seq.length())+k];
			}
		      //		      cout<<"score = "<<score<<endl;

		      int score2=score+.5;
		      //		      cout<<"score2 = "<<score2<<endl;
		      //		      corr_array->increment_element(i,j,(int)(score+.5));
		      corr_array->increment_element(i,j,score2);
		      //		      cout<<"Incremented the corr_array..."<<endl;
		      
		      if(seq.length()>=mbkl)
			{
			  bitterv[get_bitterv_index(i,j)].add_new_score(score);
			}
		      //		      cout<<"got through the if block..."<<endl;
		      
		      gh_index2++;
		    }
		}
	      gh_index1++;
	    }
	  //	  cout<<"done in the next loop"<<endl;
	}//end of new code

      //      cout<<"out of one of the loops. Does it get here?"<<endl;

      
      /*
      for(int i=0; i<number_of_genomes; i++)
	{
	  if(gh[i])
	    {
	      for(int j=i+1; j<number_of_genomes; j++)
		{
		  if(gh[j])
		    {
		      double score = get_score(i,j,seq);
		      corr_array->increment_element(i,j,(int)(score+.5));
		      
		      if(seq.length()>=mbkl)
			{
			  bitterv[get_bitterv_index(i,j)].add_new_score(score);
			} 
		    }
		}
	    }
	}

      */
      /*      set<int>::iterator it1;
      set<int>::iterator it2;      
      
      for(it1=gnum_holder.begin(); it1!=gnum_holder.end(); it1++)
	{
	  for(it2=gnum_holder.begin(); it2!=gnum_holder.end(); it2++)
	    {
	      int ig=*it1; //ig is a gnum
	      int jg=*it2; //jg is a gnum
	      if(ig<jg)
		{
		  double score = get_score(ig,jg,seq);
		  corr_array->increment_element(ig,jg,(int)(score+.5));

		  if(seq.length()>=mbkl)
		    {
		      bitterv[get_bitterv_index(ig,jg)].add_new_score(score);
		    }
		}
	    }
	}
      */
      //fill in the diagonal.
      /*      set<int>::iterator it3;
      for(it3=gnum_diag_hit.begin(); it3!=gnum_diag_hit.end(); it3++)
	{
	  int it3_int = *it3;
	  double score = get_score(it3_int, it3_int,seq);
	  corr_array->increment_element(it3_int,it3_int,(int)(score+.5));

	  }*/

      delete[] score_array;
    }

  delete[] gh;

  if(MCTR_DEBUG)
    {
      cout<<"End of increment_bits_it0(...)"<<endl;
    }
}

int mctr::set_gh(bool gh[], vector<int>& gnums, int first_row, int last_row)
{
  int count=0;
  for(int i=first_row; i<=last_row; i++)
    {
      if(!gh[gnums[i]])
	{
	  gh[gnums[i]]=true;
	  count++;
	}
    }
  return count;
}

void mctr::set_score_array(double * score_array,string seq,int count,int gen_ord)
{
  for(int j=0; j<seq.length() && seq[j]!='^'; j++)
    {
      score_array[count*seq.length()+j]=bitsv[gen_ord].d[get_aa_int[(int)seq[j]-65]];
    }
}

/*******************************************
increment_bits_it1(...)
Use when running a problematic (HGT) pair against a reference set.
run_type==1
 - The 2 problem proteomes are in positions 0 and 1 (or -2,-1). 
 - Ignore any groups of sequences in which these two proteomes are both absent.
*********************************************/
void mctr::increment_bits_it1(int first_row, int last_row,string seq)//,ofstream& alalstream, ofstream& alstream)
{
  set<int> gnum_holder; //using a set to deal with repeats within the same organism (sets are non-redundant)
  bool store_gene_ids=false; //Used for storing gene numbers in bitterv because their score >= match_score_cutoff. Even when msc is true,

  bool full=false;
  set<int> gnum_diag_hit;
  pair<set<int>::iterator,bool> ret;

  //Store all ordinal IDs in a set.
  for(int i=first_row; i<=last_row && !full; i++)
    {
      //      int gnum=gnums.at(i);
      int gnum=gnums[i];
      ret = gnum_holder.insert(gnum);
      if(gnum_holder.size()>=number_of_genomes && gnum_diag_hit.size()>=number_of_genomes)
	{
	  full=true;
	}
      if(ret.second==false) //This means gnum was already in the set. Since the sequence is present twice in the same organism, that means the organism has a match within itself; therefore add the number to the diagonal set.
	{
	  gnum_diag_hit.insert(gnum);
	}
    }

  //Case 1: No matches between organisms or within organisms.
  if(gnum_holder.size()<=1 && gnum_diag_hit.size()==0)
    {
      //do nothing
    }

  //Case 2: A match within an organism and not between organisms
  else if(gnum_holder.size()<=1 && gnum_diag_hit.size()>=1) 
    {
      set<int>::iterator it3;
      for(it3=gnum_diag_hit.begin(); it3!=gnum_diag_hit.end(); it3++)
	{
	  int it3_int = *it3;
	  double score = get_score(it3_int,it3_int,seq);
	  corr_array->increment_element(it3_int,it3_int,(int)(score+.5));
	}
    }

  //Case 3: A match between ALL organisms
  else if(gnum_holder.size()>=number_of_genomes) //shortcut
    {
      for(int i=0; i<number_of_genomes-1; i++)
	{
	  for(int j=i+1; j<number_of_genomes; j++)
	    {
	      double score = get_score(i,j,seq);
	      corr_array->increment_element(i,j,(int)(score+.5));
	      if(seq.length()>=mbkl)
		{
		  bitterv[get_bitterv_index(i,j)].add_new_score(score);
		}
	      
	      if( (i<=1 || j<=1) ) //if either problem org is represented.
		{
		  if(seq.length()>=HGT_MIN_LENGTH)
		    {
		      store_gene_ids=true;
		    }
		}
	    }
	}
      
      if(gnum_diag_hit.size()>0)
	{
	  set<int>::iterator it1;
	  for(it1=gnum_diag_hit.begin(); it1!=gnum_diag_hit.end(); it1++)
	    {
	      int it1_int = *it1;
	      double score = get_score(it1_int,it1_int,seq);
	      corr_array->increment_element(it1_int,it1_int,(int)(score+.5));
	    }
	}
    }

  //Case 4: A match between SOME non-identical organisms
  else
    {
      set<int>::iterator it1;
      set<int>::iterator it2;      
      
      for(it1=gnum_holder.begin(); it1!=gnum_holder.end(); it1++)
	{
	  for(it2=gnum_holder.begin(); it2!=gnum_holder.end(); it2++)
	    {
	      int ig=*it1; //ig is a gnum
	      int jg=*it2; //jg is a gnum
	      if(ig<jg)
		{
		  double score = get_score(ig,jg,seq);

		  corr_array->increment_element(ig,jg,(int)(score+.5));
		  if(seq.length()>=mbkl)
		    {
		      bitterv[get_bitterv_index(ig,jg)].add_new_score(score);
		    }

		  if(ig<=1 || jg<=1)
		    {
		      if(seq.length()>=HGT_MIN_LENGTH)
			{
			  store_gene_ids=true;
			}
		    }
		}
	    }
	}

      set<int>::iterator it3;
      for(it3=gnum_diag_hit.begin(); it3!=gnum_diag_hit.end(); it3++)
	{
	  int it3_int = *it3;
	  double score = get_score(it3_int, it3_int,seq);
	  corr_array->increment_element(it3_int,it3_int,(int)(score+.5));
	}
    }

  if(store_gene_ids)
    {
      set<int>::iterator it1;
      set<int>::iterator it2;
      int bottom=-1; //Don't let repeats get counted; once a gnum has been handled, ignore all the rest of the same gnum. Fortunately, the vector is sorted.
      bool hgt_match=true;

      for(int i=first_row; i<last_row; i++)
	{
	  for(int j=i+1; j<=last_row; j++)
	    {
	      int ig = gnums[i]; //ig is the NON-ORIGINAL genome ordinal.
	      int jg = gnums[j]; //jg is the next NON-ORIGNAL genome ordinal (equal or larger than ig)
	      int gig = gene_ids[i];
	      int gjg = gene_ids[j];

	      if(ig!=jg) //If the matches are not from the same genome
		{
		  double score = get_score(ig,jg,seq);
		  int index=get_bitterv_index(ig,jg);
		  
		  if(score>=bitterv[index].get_msc() || hgt_match)
		    {
		      std::pair<int,int> hit1;
		      std::pair<int,int> hit2;
		      hit1=std::make_pair(ig,gig); //proteome 1 ids
		      hit2=std::make_pair(jg,gjg); //proteome 2 ids
		      if( (ig==0 && jg==1) || (ig==1 && jg==0) ) //An entry for the HGT map.
			{
			  if(hgt_map.find(std::pair<int,int>(ig,gig))!=hgt_map.end())
			    {
			      std::map<std::pair<int,int>,int>::iterator it;
			      it = hgt_map.find(std::pair<int,int>(ig,gig));
			      (it->second)++;
			    }
			  //entry doesn't exist.
			  else
			    {
			      std::pair<std::pair<int,int>,int> foo;
			      foo = std::make_pair(hit1,1); //The 1 stands for the starting score, i.e. 1 hit.
			      hgt_map.insert(foo);			      
			    }
			  //second entry exists in the hgt map
			  if(hgt_map.find(std::pair<int,int>(jg,gjg))!=hgt_map.end())
			    {
			      std::map<std::pair<int,int>,int>::iterator it;
			      it = hgt_map.find(std::pair<int,int>(jg,gjg));
			      (it->second)++;
			    }
			  else
			    {
			      std::pair<std::pair<int,int>,int> foo;
			      foo = std::make_pair(hit2,1);
			      hgt_map.insert(foo);
			    }
			} //end of ig==0 and jg==1
		      
		      //Not a match between the problems; a match with a ref.
		      else if(ig==0 || ig==1)
			{
			  if(con_map.find(std::pair<int,int>(ig,gig))!=con_map.end())
			    {
			      std::map<std::pair<int,int>,int>::iterator it;
			      it = con_map.find(std::pair<int,int>(ig,gig));
			      (it->second)++;
			    }
			  //entry doesn't exist.
			  else
			    {
			      std::pair<std::pair<int,int>,int> foo;
			      foo = std::make_pair(hit1,1); //The 1 stands for the starting score, i.e. 1 hit.
			      con_map.insert(foo);			      
			    }
			}
		      else if(jg==0 || jg==1)
			{
			  //second entry exists in the con map
			  if(con_map.find(std::pair<int,int>(jg,gjg))!=con_map.end())
			    {
			      std::map<std::pair<int,int>,int>::iterator it;
			      it = con_map.find(std::pair<int,int>(jg,gjg));
			      (it->second)++;
			    }
			  else
			    {
			      std::pair<std::pair<int,int>,int> foo;
			      foo = std::make_pair(hit2,1);
			      con_map.insert(foo);
			    } 
			} //end of jg==0 or jg==1
		    } //end of score > msc
		} //end of ig<jg and ig>bottom
	    } //end of first for
	} //end of second for
    } //end of if(msc)
}

/*******************************************
increment_bits_it2(...)
- When running just the HGT pair, HGT-genes removed.
*********************************************/
void mctr::increment_bits_it2(int first_row, int last_row,string seq)//,ofstream& alalstream, ofstream& alstream)
{
  set<int> gnum_holder; //using a set to deal with repeats within the same organism (sets are non-redundant)
  bool store_gene_ids=false; //Used for storing gene numbers in bitterv because their score >= match_score_cutoff. Even when msc is true,

  bool full=false;
  set<int> gnum_diag_hit;
  pair<set<int>::iterator,bool> ret;

  //Store all ordinal IDs in a set.
  for(int i=first_row; i<=last_row && !full; i++)
    {
      int gnum=gnums[i];
      ret = gnum_holder.insert(gnum);
      if(gnum_holder.size()>=number_of_genomes && gnum_diag_hit.size()>=number_of_genomes)
	{
	  full=true;
	}
      if(ret.second==false) //This means gnum was already in the set. Since the sequence is present twice in the same organism, that means the organism has a match within itself; therefore add the number to the diagonal set.
	{
	  gnum_diag_hit.insert(gnum);
	}
    }

  if(MCTR_DEBUG)
    {
      cout<<"In increment_bits_it2(...)"<<endl;
      cout<<"first_row = "<<first_row<<" last_row = "<<last_row<<" seq = "<<seq<<endl;
      cout<<"ret = ";
      set<int>::iterator it;
      for(it=gnum_holder.begin(); it!=gnum_holder.end(); it++)
	{
	  cout<<*it<<" ";
	}
      cout<<endl;
    }


  //cluster_array[gnum_holder.size()][(int)(score+.5)]++;

  //Case 1: No matches between organisms or within organisms.
  if(gnum_holder.size()<=1 && gnum_diag_hit.size()==0)
    {
      //do nothing
    }
  //Case 2: A match within an organism and not between organisms
  else if(gnum_holder.size()<=1 && gnum_diag_hit.size()>=1) 
    {
      set<int>::iterator it3;
      for(it3=gnum_diag_hit.begin(); it3!=gnum_diag_hit.end(); it3++)
	{
	  int it3_int = *it3;
	  double score = get_score(it3_int,it3_int,seq);
	  corr_array->increment_element(it3_int,it3_int,(int)(score+.5));
	}
    }

  //Case 3: A match between ALL organisms
  else if(gnum_holder.size()>=number_of_genomes) //shortcut
    {
      for(int i=0; i<number_of_genomes-1; i++)
	{
	  for(int j=i+1; j<number_of_genomes; j++)
	    {
	      double score = get_score(i,j,seq);
	      corr_array->increment_element(i,j,(int)(score+.5));
	      if(seq.length()>=mbkl)
		{
		  //		  bitterv.at(get_bitterv_index(i,j)).add_new_score(score);
		  bitterv[get_bitterv_index(i,j)].add_new_score(score);
		}

	      //	      if(msc && (i<=1 || j<=1) && score>=bitterv.at(get_bitterv_index(i,j)).get_msc() && bitterv.at(get_bitterv_index(i,j)).get_msc()>-1)
	      if(msc && (i<=1 || j<=1) && score>=bitterv[get_bitterv_index(i,j)].get_msc() && bitterv[get_bitterv_index(i,j)].get_msc()>-1)
		{
		  if(seq.length()>=6)
		    {
		      store_gene_ids=true;
		    }
		}
	    }
	}
      
      if(gnum_diag_hit.size()>0)
	{
	  set<int>::iterator it1;
	  for(it1=gnum_diag_hit.begin(); it1!=gnum_diag_hit.end(); it1++)
	    {
	      int it1_int = *it1;
	      double score = get_score(it1_int,it1_int,seq);
	      corr_array->increment_element(it1_int,it1_int,(int)(score+.5));
	    }
	}
    }

  //Case 4: A match between SOME non-identical organisms
  else
    {
      set<int>::iterator it1;
      set<int>::iterator it2;      
      
      for(it1=gnum_holder.begin(); it1!=gnum_holder.end(); it1++)
	{
	  for(it2=gnum_holder.begin(); it2!=gnum_holder.end(); it2++)
	    {
	      int ig=*it1; //ig is a gnum
	      int jg=*it2; //jg is a gnum
	      if(ig<jg)
		{
		  double score = get_score(ig,jg,seq);

		  corr_array->increment_element(ig,jg,(int)(score+.5));
		  if(seq.length()>=mbkl)
		    {
		      //		      bitterv.at(get_bitterv_index(ig,jg)).add_new_score(score);
		      bitterv[get_bitterv_index(ig,jg)].add_new_score(score);
		    }
		  if(msc && (ig<=1 || jg<=1) && score>=bitterv[get_bitterv_index(ig,jg)].get_msc() && bitterv[get_bitterv_index(ig,jg)].get_msc()>-1)
		    {
		      if(seq.length()>=6)
			{
			  store_gene_ids=true;
			}
		    }
		}
	    }
	}

      set<int>::iterator it3;
      for(it3=gnum_diag_hit.begin(); it3!=gnum_diag_hit.end(); it3++)
	{
	  int it3_int = *it3;
	  double score = get_score(it3_int, it3_int,seq);
	  corr_array->increment_element(it3_int,it3_int,(int)(score+.5));
	}
    }
} // end of increment_bits_it2(...)

int mctr::get_aa(char aa)
{
  if(aa=='A')
    return 0;
  else if(aa=='C')
    return 1;
  else if(aa=='D')
    return 2;
  else if(aa=='E')
    return 3;
  else if(aa=='F')
    return 4;
  else if(aa=='G')
    return 5;
  else if(aa=='H')
    return 6;
  else if(aa=='I')
    return 7;
  else if(aa=='K')
    return 8;
  else if(aa=='L')
    return 9;
  else if(aa=='M')
    return 10;
  else if(aa=='N')
    return 11;
  else if(aa=='P')
    return 12;
  else if(aa=='Q')
    return 13;
  else if(aa=='R')
    return 14;
  else if(aa=='S')
    return 15;
  else if(aa=='T')
    return 16;
  else if(aa=='V')
    return 17;
  else if(aa=='W')
    return 18;
  else if(aa=='Y')
    return 19;
  else if(aa=='^')
    return 20;
  else
    {
      cout<<"Invalid character: "<<aa<<endl;
      exit(1);
      return -1;
    }
}

char mctr::get_aa_char_fn(int i)
{
  if(i<0 || i>=20)
    {
      cout<<"Bad value for int i in fn get_aa_char_fn(int): i = "<<i<<endl<<"Exiting."<<endl;
      exit(1);
    }
  else if(i==0) return 'A';
  else if(i==1) return 'C';
  else if(i==2) return 'D';
  else if(i==3) return 'E';
  else if(i==4) return 'F';
  else if(i==5) return 'G';
  else if(i==6) return 'H';
  else if(i==7) return 'I';
  else if(i==8) return 'K';
  else if(i==9) return 'L';
  else if(i==10) return 'M';
  else if(i==11) return 'N';
  else if(i==12) return 'P';
  else if(i==13) return 'Q';
  else if(i==14) return 'R';
  else if(i==15) return 'S';
  else if(i==16) return 'T';
  else if(i==17) return 'V';
  else if(i==18) return 'W';
  else if(i==19) return 'Y';
  else
    {
      cout<<"Invalid value again in get_aa_char_fn(int)."<<endl;
      exit(1);
    }
    
}

/******************************
Write out data for a normal, main run.
 *****************************/
void mctr::write_out_data_it0()
{
  cout<<"mctr::write_out_data_it0()"<<endl;

  //print values to screen.
  for(int i=0; i<bitterv.size(); i++)
    {
      //      bitterv.at(i).write_out_values();
      bitterv[i].write_out_values();
    }

  //Write all plots to data.txt in DATA/
  ofstream dstream;
  string datafile = path+"/"+bfn+"/"+data_directory+"/data.txt";
  u.open_ofile_ne(dstream,datafile,"datafile");
  dstream<<max_bit_sum<<endl<<endl;
  
  for(int i=0; i<number_of_genomes-1; i++)
    {
      for(int j=i+1; j<number_of_genomes; j++)
	{
	  dstream<<"ids "<<i<<" "<<j<<endl;
	  //dstream<<"mean "<<bitterv.at(get_bitterv_index(i,j)).get_average()<<endl;
	  //dstream<<"rms "<<bitterv.at(get_bitterv_index(i,j)).get_rms()<<endl;
	  //dstream<<"min_bitscore "<<bitterv.at(get_bitterv_index(i,j)).get_min()<<endl;	  

	  dstream<<"mean "<<bitterv[get_bitterv_index(i,j)].get_average()<<endl;
	  dstream<<"rms "<<bitterv[get_bitterv_index(i,j)].get_rms()<<endl;
	  dstream<<"min_bitscore "<<bitterv[get_bitterv_index(i,j)].get_min()<<endl;

	  int holder_size = max_bit_sum;
	  int holder[holder_size];
	  for(int k=0; k<holder_size; k++)
	    {
	      holder[k]=corr_array->get_element(i,j,k);	      
	    }
	  for(int k=0; k<holder_size; k++)
	    {
	      dstream<<holder[k]<<endl;
	    }
	  dstream<<endl;
	}
      dstream<<endl;
    }
  dstream.close(); //Done writing out plots.
  
  ofstream diagstream;
  string outfile4=path+"/"+bfn+"/"+data_directory+"/diagonals.txt";
  u.open_ofile_ne(diagstream,outfile4,"outfile4");
  
  for(int i=0; i<number_of_genomes; i++)
    {
      diagstream<<i<<" to "<<i<<endl;
      for(int j=0; j<max_bit_sum; j++)
	{
	  diagstream<<corr_array->get_element(i,i,j)<<endl;
	}
      diagstream<<endl;
    }
  diagstream.close();
} //end of write_out_data_it0()

/******************************
Write out data for run between problematic (HGT) pair and reference set.
- This does not include writing out plots.
- This only includes writing out genes that are suspicious for each pair.
 *****************************/
void mctr::write_out_data_it1()
{
  genes_found=0;

  ofstream con_stream;
  string con_outfile = path+"/"+bfn+"/con_msc_for_filtering.txt";
  u.open_ofile_app(con_stream,con_outfile,"con_outfile");
  
  ofstream hgt_stream;
  string hgt_outfile = path+"/"+bfn+"/hgt_msc_for_filtering.txt";
  u.open_ofile_app(hgt_stream,hgt_outfile,"hgt_outfile");
  
  con_stream<<"Ordinals: "<<prob1<<" "<<prob2<<endl;
  hgt_stream<<"Ordinals: "<<prob1<<" "<<prob2<<endl;
  
  for(std::map<std::pair<int,int>, int>::iterator it=con_map.begin(); it!=con_map.end(); ++it)
    {
      con_stream<<(it->first).first<<" "<<(it->first).second<<" "<<it->second<<endl;
    }
  
  for(std::map<std::pair<int,int>, int>::iterator it=hgt_map.begin(); it!=hgt_map.end(); ++it)
    {
      hgt_stream<<(it->first).first<<" "<<(it->first).second<<" "<<it->second<<endl;
    }
  
  //Go through all the genes that shared long matches between the two problem organisms.
  hgt_stream<<endl<<"HGT_instances: "<<endl;
  for(std::map<std::pair<int,int>, int>::iterator hgt_it=hgt_map.begin(); hgt_it!=hgt_map.end(); ++hgt_it)
    {
      bool in_there=false;
      int in_there_count=0;
      for(std::map<std::pair<int,int>, int>::iterator it=con_map.begin(); it!=con_map.end(); ++it)
	{
	  if(hgt_it->first.first==it->first.first && hgt_it->first.second==it->first.second)
	    {
	      in_there=true;
	      in_there_count=it->second;
	    }
	}

      int ref_hits=0; 
      int ref_hits_other=0;

      if(!in_there || in_there_count<=CON_MAP_CUTOFF)
	{
	  cout<<"/nMismatch between conserved set and hgt pair for "<<hgt_it->first.first<<" "<<hgt_it->first.second<<" "<<hgt_it->second<<endl;
	  cout<<"in_there="<<in_there<<endl;
	  cout<<"in_there_count="<<in_there_count<<endl;
	  genes_found++;
	  hgt_stream<<hgt_it->first.first<<" "<<hgt_it->first.second<<" "<<hgt_it->second<<endl;
	  if(hgt_it->first.first==0)
	    {
	      cout<<"Info line = "<<info_lines0[hgt_it->first.second]<<endl;
	      hgt_stream<<"Info line = "<<info_lines0[hgt_it->first.second]<<endl<<proteins0[hgt_it->first.second]<<endl;
	    }
	  else if(hgt_it->first.first==1)
	    {
	      cout<<"Info line = "<<info_lines1[hgt_it->first.second]<<endl;
	      hgt_stream<<"Info line = "<<info_lines1[hgt_it->first.second]<<endl<<proteins1[hgt_it->first.second]<<endl;
	    }
	}
    }  
  
  hgt_stream<<"End_of_hgt_genes"<<endl<<endl;      
  con_stream.close();
  hgt_stream.close();
} //end of write_out_data_it1()


/******************************
When running just a pair, HGT-genes removed in previous run against ref set.
- Write out new plot
 *****************************/
void mctr::write_out_data_it2()
{
  if(MCTR_DEBUG)
    {
      cout<<"void mctr::write_out_data_it2()"<<endl;
    }

#pragma omp critical
  if(run_type==2)
    {
      cout<<"Sort of writing it out... Here it is for now, with a staller at the end..."<<endl;
      cout<<"Writing out the data...."<<endl;
      ofstream outstream2;
      outstream2.open( (path+"/"+bfn+"/"+data_directory+"/data_after_hgt_cleanup.txt").c_str(), ios::out | ios::app);
      while(outstream2.fail())
	{
	  cout<<"Failed to open outfile. Does directory DATA exist?"<<endl;
	  int staller;
	  cout<<"Stalling. Create the damn directory please: (any int) "<<endl;
	  cin>>staller;
	  outstream2.open( (path+"/"+bfn+"/"+data_directory+"/data_after_hgt_cleanup.txt").c_str());
	}
      cout<<"Successfully opened outfile to data.txt"<<endl;

      outstream2<<max_bit_sum<<endl<<endl;

      for(int i=0; i<number_of_genomes-1; i++)
	{
	  for(int j=i+1; j<number_of_genomes; j++)
	    {
	      outstream2<<prob1<<" to "<<prob2<<endl;

	      //	      cout<<bitterv.at(get_bitterv_index(i,j)).get_average()<<endl;
	      //	      outstream2<<"mean = "<<bitterv.at(get_bitterv_index(i,j)).get_average()<<endl;
	      //	      outstream2<<"rms = "<<bitterv.at(get_bitterv_index(i,j)).get_rms()<<endl;
	      //	      outstream2<<"min_bitscore = "<<bitterv.at(get_bitterv_index(i,j)).get_min()<<endl;

	      cout<<bitterv[get_bitterv_index(i,j)].get_average()<<endl;
	      outstream2<<"mean = "<<bitterv[get_bitterv_index(i,j)].get_average()<<endl;
	      outstream2<<"rms = "<<bitterv[get_bitterv_index(i,j)].get_rms()<<endl;
	      outstream2<<"min_bitscore = "<<bitterv[get_bitterv_index(i,j)].get_min()<<endl;

	      int holder_size = max_bit_sum;
	      int holder[holder_size];
	      for(int k=0; k<holder_size; k++)
		{
		  holder[k]=corr_array->get_element(i,j,k);	      
		}
	      for(int k=0; k<holder_size; k++)
		{
		  outstream2<<holder[k]<<endl;
		}
	      outstream2<<endl;
	    }
	  outstream2<<endl;
	}
      outstream2.close();      
    }

  cout<<"Bottom of write out pair run data"<<endl;
}
  
/*************************

 ************************/
void mctr::read_in_files()
{
  if(MCTR_DEBUG)
    {
      cout<<"mctr::read_in_files()"<<endl;
    }

  DIR *dp;
  struct dirent *dirp;
  string pathfile=path+"/"+bfn+"/"+merged_directory+"/";
  cout<<"pathfile = "<<pathfile<<endl;

  if((dp = opendir( (pathfile.c_str()))) == NULL)
    {
      cout<<"Failed to open path to "<<pathfile<<endl;
      exit(1);
    }
  
  while ((dirp = readdir(dp)) != NULL) 
    {
      string temp = string(dirp->d_name);
      if(temp.length()>5)
	{
	  files.push_back(string(dirp->d_name));
	}
      else
	{
	  cout<<"Rejecting file named "<<temp<<endl;
	}
    }

  closedir(dp);
  if(MCTR_DEBUG)
    {
      cout<<"mctr::mctr() : files.size() = "<<files.size()<<endl;
      for(int i=0; i<files.size(); i++)
	{
	  //	  cout<<files.at(i)<<endl;
	  cout<<files[i]<<endl;
	}
    }
}

int mctr::calculate_max_bit_sum()
{
  return 2*max_single_bit_value*tag_length; //2* because the bit values are halved...
}

void mctr::open_info_file(ifstream& infostream)
{
  string infile=path+"/"+bfn+"/"+bfn+"_info.txt";
  u.open_ifile(infostream,infile,"infile");
  int nog_ref,nog,nog_total;
  infostream>>nog_ref>>nog>>tag_length;
  number_of_genomes=nog+nog_ref;
  number_of_genomes_ref=nog_ref; 
}

/************************************************************************************************
For a second or third iteration run, when a pair has been selected that may haveHGT and is either being run against a reference set or just run against each other with the suspicious genes marked for exclusion, THIS FUNCTION READS THE INFO FILE AND ADDS ALL .FAA FILES ASSOCIATED WITH EACH OF THE PROBLEMATIC ORGANISMS ONTO A PAIR OF STRING VECTORS.
 ************************************************************************************************/
void mctr::read_in_probfiles(ifstream& infostream,vector<string>& files1, vector<string>& files2)
{
  if(MCTR_DEBUG)
    {
      cout<<"In mctr::read_in_probfiles(...)"<<endl;
      cout<<"files1.size() = "<<files1.size()<<endl;
      cout<<"files2.size() = "<<files2.size()<<endl;
      cout<<"prob1="<<prob1<<" prob2="<<prob2<<endl;
    }

  int ord=-1;
  string eater;

  infostream>>eater;
  infostream>>eater;
  infostream>>eater;
  infostream>>ord;  

  if(MCTR_DEBUG)
    {
      cout<<"In mctr::read_in_probfiles(...): read in ord="<<ord<<endl;
    }

  while(!infostream.eof())
    {
      string org_name;
      string full_path;
      infostream>>org_name>>full_path;
      if(ord==prob1)
	{
	  files1.push_back(full_path);
	}
      else if(ord==prob2)
	{
	  files2.push_back(full_path);
	}
      infostream>>ord;
      if(MCTR_DEBUG)
	{
	  cout<<"In mctr::read_in_probfiles(...): read in ord="<<ord<<endl;
	}
    }
}

//For the 2nd iteration, in which a pair flagged for HGT is run against a reference set to identify genes with suspiciously long kmer matches.
// - infostream -> <bfn>_info.txt [the original info file]
void mctr::setup_it1(ifstream& infostream)
{
  vector<string> files1;
  vector<string> files2;      
  read_in_probfiles(infostream,files1,files2); //just get the faa files for each problem organism.

  get_problem_tags(files1,"0");
  get_problem_tags(files2,"1");


    std::sort(problem_tags.begin(), problem_tags.end());


  if(MCTR_DEBUG)
    {
      cout<<"void mctr::setup_it1(ifstream& infostream)"<<endl;
      cout<<"problem_tags.size()="<<problem_tags.size()<<endl;
      for(int i=0; i<10 && i<problem_tags.size(); i++)
	{
	  //	  cout<<prob1<<" "<<prob2<<" "<<problem_tags.at(i)<<endl;
	  cout<<prob1<<" "<<prob2<<" "<<problem_tags[i]<<endl;
	}
    }
}

/***************************************
- Read in the tags for problem proteomes 0 and 1
- 
 **************************************/
void mctr::setup_it2(ifstream& infostream)
{
  vector<string> files1;
  vector<string> files2;      
  read_in_probfiles(infostream,files1,files2);

  set<int> hgt_genes0;
  set<int> hgt_genes1;

  string eater;

  ifstream instream;
  instream.open( (path+"/"+bfn+"/hgt_msc_for_filtering.txt").c_str());
  if(instream.fail())
    {
      cout<<"Failed to open hgt_msc_for_filtering.txt"<<endl;
      exit(1);

    }
  else
    {
      cout<<"Successfully opened hgt_msc_for_filtering.txt"<<endl;
    }
  
  int e1=-1;
  int e2=-1;
  
  eater="";
  
  while( (e1!=prob1 || e2!=prob2) && !instream.eof() )
    {
      while(eater!="Ordinals:" && !instream.eof() )
	{
	  instream>>eater;
	}
      instream>>e1>>e2;
      eater="";
    }
  
  while(eater!="HGT_instances:" && !instream.eof() )
    {
      instream>>eater;
    }
  //at this point, instream should be poised at hgt genes.
  
  instream>>eater; //should eat a 0 or 1 OR "End_of_hgt_genes" if there are no hgt genes!!!
  
  while(eater!="End_of_hgt_genes")
    {
      int gene_id,count;
      instream>>gene_id>>count;
      char buffer[2000];
      string info_line,p;
      instream>>info_line;
      instream.getline(buffer,5000);
      info_line+=buffer;
      instream>>p;
      if(p=="" || info_line=="")
	{
	  cout<<"p is empty in mctr."<<endl;
	  exit(1);
	}
      
      if(eater=="0")
	{
	  hgt_genes0.insert(gene_id);
	}
      else if(eater=="1")
	{
	  hgt_genes1.insert(gene_id);
	}
      instream>>eater;
    }
  
  
  instream.close();
    
  get_problem_tags(files1,"0",hgt_genes0);
  get_problem_tags(files2,"1",hgt_genes1);
  std::sort(problem_tags.begin(), problem_tags.end());
}

void mctr::create_directories()
{
  int system_return = system ( ("rm -r "+path+"/"+bfn+"/"+data_directory).c_str());
  system_return = system ( ("mkdir "+path+"/"+bfn+"/"+data_directory).c_str());
}

double mctr::get_score(int i1,int i2, string seq)
{
  double result=0;
  for(int i=0; i<seq.length(); i++)
    {
      //      int index = get_aa(seq[i]);
      int index = get_aa_int[(int)seq[i]-65];

      if(index>=0 && index<20)
	{
	  //	  result=result+bitsv.at(i1).d[get_aa(seq[i])]+bitsv.at(i2).d[get_aa(seq[i])];
	  //result=result+bitsv.at(i1).d[get_aa_int[(int)seq[i]-65]]+bitsv.at(i2).d[get_aa_int[(int)seq[i]-65]];
	  result=result+bitsv[i1].d[get_aa_int[(int)seq[i]-65]]+bitsv[i2].d[get_aa_int[(int)seq[i]-65]];
	}
      else
	{
	  cout<<"index out of bounds in get_score(...). seq = "<<seq<<endl;
	  exit(1);
	}

    }
  return result;

}

void mctr::populate_bitterv()
{
  if(run_type==1) //pair against ref_set
    {
      int total=refsetsize;
      if(prob1>=refsetsize)
	total++;
      if(prob2>=refsetsize)
	total++;

      for(int i=0; i<total-1; i++)
	{
	  for(int j=i+1; j<total; j++)
	    {
	      bitter btemp;
	      btemp.set_ids(i,j);
	      btemp.set_msc(-1);
	      bitterv.push_back(btemp);
	    }
	}
    }
  else if(run_type==2) //pair against itself only
    {
      for(int i=0; i<number_of_genomes-1; i++)
	{
	  for(int j=i+1; j<number_of_genomes; j++)
	    {
	      bitter btemp;
	      btemp.set_ids(i,j);
	      btemp.set_msc(-1);
	      bitterv.push_back(btemp);
	    }
	}
    }
  else //normal case.
    {
      for(int i=0; i<number_of_genomes-1; i++)
	{
	  for(int j=i+1; j<number_of_genomes; j++)
	    {
	      bitter btemp;
	      btemp.set_ids(i,j);
	      btemp.set_msc(-1);
	      bitterv.push_back(btemp);
	    }
	}
    }
}

int mctr::get_bitterv_index(int i1, int i2)
{
  if(run_type==0)
    {
      if(i1<i2)
	{
	  int index = ((i1*number_of_genomes)+i2) - (( (i1+1)*((i1+1)+1)))/2;
	  return index;
	}
      else if(i2<i1)
	{
	  int index = ((i2*number_of_genomes)+i1) - (((i2+1)*((i2+1)+1)))/2;
	  return index;
	}
      else
	{
	  cout<<"i1 = "<<i1<<" i2 = "<<i2<<endl;
	  cout<<"i1 = i2 in mctr::get_bitterv_index(...)"<<endl;
	  exit(1);
	}      
    }
  else 
    {
      if(i1<i2)
	{
	  int index = ((i1*number_of_genomes_ref)+i2) - (( (i1+1)*((i1+1)+1)))/2;
	  //cout<<"i1 = "<<i1<<" i2 = "<<i2<<" index = "<<index<<endl;
	  return index;
	}
      else if(i2<i1)
	{
	  int index = ((i2*number_of_genomes_ref)+i1) - (((i2+1)*((i2+1)+1)))/2;
	  return index;
	}
      else
	{
	  cout<<"i1 = "<<i1<<" i2 = "<<i2<<endl;
	  cout<<"i1 = i2 in mctr::get_bitterv_index(...)"<<endl;
	  exit(1);
	}
    }
}

// For iteration 1.
void mctr::get_problem_tags(vector<string> files,string ordinal)
{
  int gene_id=0;

  for(int i=0; i<files.size(); i++)
    {
      ifstream instream;
      //      instream.open( files.at(i).c_str());
      instream.open( files[i].c_str());
      if(instream.fail())
	{
	  //	  cout<<"Failed to open "<<files.at(i)<<" in mctr::get_problem_tags(vector<string>)"<<endl;
	  cout<<"Failed to open "<<files[i]<<" in mctr::get_problem_tags(vector<string>)"<<endl;
	  exit(1);
	}
      
      string eater;
      instream>>eater;
      while(!instream.eof())
	{
	  if(eater[0]=='>')
	    {
	      char buffer[2000];
	      instream.getline(buffer,2000);
	      string info_line=eater+" "+buffer;
	      if(ordinal=="0")
		{
		  info_lines0.push_back(info_line);		  
		}
	      else
		{
		  info_lines1.push_back(info_line);
		}
	      instream>>eater;
	    }
	  string p;
	  while(!instream.eof() && eater[0]!='>')
	    {
	      p+=eater;
	      instream>>eater;
	    }

	  //at this point we have the protein.
	  if(ordinal=="0")
	    {
	      proteins0.push_back(p);
	    }
	  else
	    {
	      proteins1.push_back(p);
	    }

	  for(int q=0; q<p.length(); q++)
	    {
	      string t = p.substr(q,tag_length);
	      if(check_composition(t) && t.length()==tag_length)
		{
		  stringstream ss;
		  ss<<gene_id;
		  string s = ss.str();
		  t+=" "+ordinal+" "+s;
		  problem_tags.push_back(t);
		}
	    }
	  gene_id++;
	}
      instream.close();
    }
}

void mctr::get_problem_tags(vector<string> files,string ordinal, set<int> hgt_genes)
{
  cout<<"In get_problem_tags for pair_run."<<endl;

  int gene_id=0;
  for(int i=0; i<files.size(); i++)
    {
      ifstream instream;
      //      instream.open( files.at(i).c_str());
      instream.open( files[i].c_str());
      if(instream.fail())
	{
	  //	  cout<<"Failed to open "<<files.at(i)<<" in mctr::get_problem_tags(vector<string>)"<<endl;
	  cout<<"Failed to open "<<files[i]<<" in mctr::get_problem_tags(vector<string>)"<<endl;
	  exit(1);
	}
      
      string eater;
      instream>>eater;
      while(!instream.eof())
	{
	  if(eater[0]=='>')
	    {
	      char buffer[2000];
	      instream.getline(buffer,2000);
	      string info_line=eater+" "+buffer;
	      if(ordinal=="0")
		{
		  info_lines0.push_back(info_line);		  
		}
	      else
		{
		  info_lines1.push_back(info_line);
		}
	      instream>>eater;
	    }
	  string p;
	  while(!instream.eof() && eater[0]!='>')
	    {
	      p+=eater;
	      instream>>eater;
	    }
	  //at this point we have the protein.
	  if(ordinal=="0")
	    {
	      proteins0.push_back(p);
	    }
	  else
	    {
	      proteins1.push_back(p);
	    }

	  set<int>::iterator it;
	  it = hgt_genes.find(gene_id);

	  if(it==hgt_genes.end())
	    {
	      for(int q=0; q<p.length(); q++)
		{
		  string t = p.substr(q,tag_length);
		  if(check_composition(t))
		    {
		      while(t.length()<tag_length)
			{
			  t+="^";
			}
		      stringstream ss;
		      ss<<gene_id;
		      string s = ss.str();
		      t+=" "+ordinal+" "+s;
		      problem_tags.push_back(t);
		    }
		  
		}
	    }
	  else
	    {
	    }
	  gene_id++;
	}
      instream.close();
    }
}

bool mctr::check_composition(string s)
{
  vector<int> array;
  for(int i=0; i<20; i++) 
    {
      array.push_back(0);
    }

  for(int i=0; i<s.length(); i++)
    {
      //      int index=aas.get_aa(s[i]);
      int index=get_aa_int[(int)s[i]-65];
      if(index>=0 && index<20)
	{
	  array[index]++;
	}
      else
	{
	  return false;
	}
    }

  int result=0;
  for(int i=0; i<array.size(); i++)
    {
      result+=array[i]*array[i];
    }

  if(result>sos_cutoff)
    {
      return false;
    }
  
  return true;
}

int mctr::get_sos_cutoff()
{
  return 6.5*tag_length;
}

char mctr::get_pr_lead()
{
  ifstream instream;
  string infile = path+"/"+bfn+"/TAGS/TAGS0.txt";
  instream.open(infile.c_str());
  if(instream.fail())
    {
      cout<<"In function tmerg::get_pr_lead(): Failed to open infile to "<<infile<<endl;
      exit(1);
    }
  int num_tags;
  instream>>num_tags;
  string first_tag;
  instream>>first_tag;
  return first_tag[0];
  instream.close();
}

void mctr::set_get_aa_arrays()
{
  get_aa_int[0]=0;  //A
  get_aa_int[1]=-1; //B
  get_aa_int[2]=1;  //C
  get_aa_int[3]=2;  //D
  get_aa_int[4]=3;  //E
  get_aa_int[5]=4;  //F
  get_aa_int[6]=5;  //G
  get_aa_int[7]=6;  //H
  get_aa_int[8]=7;  //I
  get_aa_int[9]=-1; //J
  get_aa_int[10]=8; //K
  get_aa_int[11]=9; //L
  get_aa_int[12]=10;//M
  get_aa_int[13]=11;//N
  get_aa_int[14]=-1;//O
  get_aa_int[15]=12;//P
  get_aa_int[16]=13;//Q
  get_aa_int[17]=14;//R
  get_aa_int[18]=15;//S
  get_aa_int[19]=16;//T
  get_aa_int[20]=-1;//U
  get_aa_int[21]=17;//V
  get_aa_int[22]=18;//W
  get_aa_int[23]=-1;//X
  get_aa_int[24]=19;//Y
  get_aa_int[25]=-1;//blank
  get_aa_int[26]=-1;//blank
  get_aa_int[27]=-1;//blank
  get_aa_int[28]=-1;//blank
  get_aa_int[29]=-1;//^

  get_aa_char[0]='A';
  get_aa_char[1]='C';
  get_aa_char[2]='D';
  get_aa_char[3]='E';
  get_aa_char[4]='F';
  get_aa_char[5]='G';
  get_aa_char[6]='H';
  get_aa_char[7]='I';
  get_aa_char[8]='K';
  get_aa_char[9]='L';
  get_aa_char[10]='M';
  get_aa_char[11]='N';
  get_aa_char[12]='P';
  get_aa_char[13]='Q';
  get_aa_char[14]='R';
  get_aa_char[15]='S';
  get_aa_char[16]='T';
  get_aa_char[17]='V';
  get_aa_char[18]='W';
  get_aa_char[19]='Y';
}